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ABSTRACT 

This  report  is  the  last  of  our  contract  Grant  FA8655-06-C-4004.  The  report  is  divided  in  three  parts. 
The  first  part  presents  a  parametric  study  of  surface  dielectric  barrier  discharges  driven  by 
nanosecond  pulses.  The  second  discusses  the  accuracy  of  the  simulation  by  comparing  two  different 
numerical  methods  used  in  this  work.  The  third  is  a  summary  of  the  results  and  conclusions  drawn 
from  the  work  performed  in  the  frame  of  this  contract. 

Attached  to  this  report  are  the  papers  related  to  this  work  and  published  during  the  course  of  the 
contract. 
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I.  SURFACE  DIELECTRIC  BARRIER  DISCHARGES  (SDBD)  DRIVEN  BY 
NANOSECOND  PULSES  :  A  PARAMETRIC  STUDY 


The  purpose  of  this  study  is  to  perform  a  parametric  study  of  the  gas  dynamics  induced  by  a 
nanosecond  surface  dielectric  barrier  discharge.  The  geometry  is  based  on  some  recent  experiments 
of  Starikovskii  and  coworkers,  [1]  (see  Fig.  1)  and  which  were  also  used  in  our  previous  report  and 
in  Ref.  [2],  Starikovskii  et  al.  introduced  a  new  type  of  surface  dielectric  barrier  discharge  using  a 
nanosecond  voltage  pulse  generator.  The  voltage  pulse  can  be  several  10s  of  kV  with  rise  and  decay 
times  on  the  order  or  less  than  10  ns.  Under  these  conditions,  the  corona  regime  that  is  responsible 
for  the  ion  wind  in  sinusoidal  regimes  is  no  longer  present  and  the  main  discharge  regime  is  a 
streamer  regime.  This  was  confirmed  by  the  quasi-zero  ion  wind  measuredin  [1],  The  authors 
showed  that  this  kind  of  discharge  was  able  to  affect  the  aerodynamic  properties  of  a  flow  along  the 
surface  in  spite  of  the  quasi-zero  ion  wind,  and  that  a  detached  flow  could  be  reattached  when  a 
nanosecond  voltage  pulse  was  applied  between  the  electrodes  at  a  repetition  rate  of  a  few  kHz  (for 
spanwise  as  well  as  streamwise  configurations  of  the  DBD  actuators  with  respect  to  the  flow 
direction).  The  computation  domain  is  recalled  in  Fig.  1. 
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Figure  1  :  Computational  domain.  The  electrodes  are  embedded  in  a  300  pm  wide  PMMA  layer 
(sr=5).  The  electrodes  are  5  mm  long  and  37  pm  wide. 

Both  the  Navier-Stokes  equations  and  plasma  equations  are  solved  on  a  few  centimeters  large 
domain  using  the  new  numerical  technique  called  Asynchronous  Adaptive  Mesh  Refinement  [4], 
[5].  The  sensitivity  to  different  parameters  such  as  maximum  voltage,  pulse  duration,  seed  electron 
density,  secondary  emission  are  presented  in  the  following  pages.  Some  prospective  results  on  the 
influence  of  an  external  airflow  on  the  discharge  are  also  given  to  point  out  tendencies. 


1.1  SENSITIVITY  TO  THE  MAXIMUM  VOLTAGE  OF  THE  NANOSECOND  PULSES 

The  most  natural  question  that  arises  when  trying  to  optimize  the  effect  on  the  neutral  gas  is 
certainly  :  does  the  effect  increase  with  voltage  ? 

This  has  been  investigated  by  performing  simulations  with  the  same  voltage  slope  during  the  rise 
and  fall  phases  of  the  pulse,  i.e.  2000  V/ns  in  the  rise  phase  and  -933  V/ns,  but  with  different  values 
of  the  voltage  maximum.  The  voltage  values  in  the  plateau  phases  are  equal  to  14,  21,  28,  35  and  42 
kV  respectively  (see  Fig.  2). 
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Figure  2  :  Shape  of  the  applied  voltage  pulses  in  the  simulations,  with  maxima  (plateau)  at  14,  21, 
28,  35  and  42  kV 

The  calculated  discharge  currents  for  the  voltage  pulses  of  Fig.  2  are  shown  in  Fig.  3.  We  see 
that  the  current  response  of  the  discharge  follows  the  same  curve  during  the  potential  rise  phase,  as 
expected.  For  28kV  and  higher  the  maximum  value  of  the  current  saturates  around  1800  A/m,  due 
to  the  fact  that  the  streamer  head  has  reached  the  maximum  extension  of  the  lower  electrode  (5 
mm).  After  that  point  there  is  still  some  heating  as  shown  in  Figs.  4a  and  5  below.  Nevertheless  this 
late  heating  could  be  optimized  by  using  a  longer  lower  electrode.  The  same  phenomenon  occurs 
during  the  second  (negative)  pulse,  for  28  kV  and  higher  voltages,  the  maximum  current  value 
appears  to  be  limited  by  the  extension  of  the  lower  electrode.  Figure  5  shows  the  total  energy 
deposition  in  the  neutral  gas  by  the  discharge.  This  deposition  seems  to  vary  quite  linearly  with  the 
voltage  amplitude.  However  the  maximum  temperature  reached  is  saturating  around  35  kV.  This 
can  be  explained  if  one  takes  into  account  that : 

-  the  maximum  temperature  is  reached  at  the  very  tip  of  the  upper  electrode 

-  the  maximum  heating  occurs  during  the  negative  current  pulse,  when  the  positive  ions  are 
collected  by  the  electrode  tip 

-  the  negative  current  pulse  occurs  around  60  ns  and  the  smallest  mesh  size  is  about  20  pm 
which  correspond  roughly  to  the  distance  covered  by  an  acoustic  wave  during  60  ns.  This  means 
that  some  energy  has  already  be  dissipated  by  the  shock  wave. 

The  maximum  temperature  could  be  further  increased  with  voltage  higher  than  35  kV  however 
this  is  only  possible  if  the  potential  variation  rate  dV/dt  is  larger. 
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Figure  3:  Calculated  discharge  currents  for  the  voltage  pulses  of  Fig.  2  ( the  maximum  voltage  is 
indicated  on  the  figure). 


Figure  4a:  Maximum  air  temperature  calculated  in  the  conditions  of  Figs.  2,  3.  The  maximum  value 
is  reached  in  a  very  small  volume  (and  is  much  smaller  when  averaged  over  a  finite  volume 
element,  e.g.  20  pm  X  20  pm  or  more,  see  Fig.  4b) 
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Figure  4b:  Maximum  air  temperature  in  the  reference  case  (14  kV)  with  different  integration 
volumes 

The  temperature  reached  by  the  hottest  cell  of  the  mesh  (about  20  pm2)  as  shown  on  Fig.  4a  is 
almost  900  K  for  a  few  ns.  However  this  data  is  difficult  to  obtain  experimentally  because  it  needs 
very  high  resolution  measurements  in  space  and  time.  The  reason  why  the  temperature  is  so  hot  in 
the  calculation  is  because  the  maximum  is  located  at  the  very  tip  of  the  electrode.  When  the 
temperature  is  averaged  over  a  greater  volume  (100  or  200  pm2  see  Fig.  4b)  the  temperature 
increase  is  on  the  order  of  100  K  which  is  comparable  to  the  measurements  of  Starikovskii  and  al. 
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Figure  5:  Total  energy  deposition  in  the  neutral  gas 
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The  voltage  amplitude  seems  to  be  a  good  design  parameter  to  obtain  some  specific  heating 
effect  on  the  flow.  The  dV/dt  is  also  an  important  parameter  that  needs  to  be  studied. 

1.2  SENSITIVITY  TO  THE  DURATION  OF  THE  NANOSECOND  PULSES 

In  the  reference  case  (14  kV  in  the  previous  section,  see  Fig.  3),  the  discharge  current  is  not  zero 
when  it  changes  sign  abruptly  shortly  after  the  voltage  starts  to  decrease,  after  the  plateau  (t=20  ns). 
We  looked  at  the  effect  of  the  duration  of  the  voltage  plateau  on  the  energy  transfer.  Pulse  durations 
from  35  to  55  ns  were  tested  as  shown  in  Fig.  6.  Indeed  it  is  possible  to  extract  more  current  during 
the  positive  pulse  (see  Fig.  7)  which  leads  to  a  higher  heating  of  the  gas  (see  Fig.  8).  Because  the 
positive  pulse  has  transferred  more  charge  to  the  surface  for  longer  pulses,  the  following  negative 
current  pulse  is  a  little  stronger.  This  explain  the  higher  temperature  obtained  (see  Fig.  9).  However 
if  the  pulse  duration  is  too  long  the  maximum  temperature  decreases  because  the  acoustic  wave  has 
started  to  dissipate  energy. 


Voltage  shape  of  the  different  pulses 


Figure  6:  Shape  of  the  applied  pulses 
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Current  response  of  the  different  pulses 


Figure  7:  Current  on  the  lower  electrode  for  different  voltage  plateau  durations  (Fig.  6) 
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Figure  8:  Total  energy  deposition  in  the  neutral  gas  for  different  voltage  plateau  durations  (Fig.  6) 
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Figure  9:  maximum  air  temperature  reached  for  different  voltage  plateau  durations  (Fig.  6) 

1.3  SENSITIVITY  TO  THE  SEED  ELECTRON  DENSITY 

Our  plasma  model  describes  the  charged  particles  with  fluid  equations.  The  discharge  is  initiated 
by  assuming  a  given,  small  and  uniform  density  of  electrons  and  ions  in  the  simulation  domain  at 
t=0.  The  effect  of  this  initial  density  is  similar  to  photoionization  since  it  provides  the  seeds 
electrons  necessary  to  ensure  the  cathode  streamer  propagation  if  the  contribution  of  secondary 
emission  from  the  surface  is  not  sufficient  or  not  fast  enough.  We  find  that  the  calculations  are  not 
very  sensitive  to  the  seed  electron  density  level.  In  any  case  the  impact  is  limited  to  the  first  half  of 
the  positive  current  pulse  (see  Figs.  10a  &  b)  and  leads  to  no  significant  variation  in  the  total  energy 
deposition  (see  Fig.  11).  We  also  performed  calculations  with  a  non  uniform  initial  density  of  seed 
electrons  (1010  m~3  in  the  near  vicinity  of  the  exposed  electrode  and  zero  everywhere  else,  see  Fig. 
10a).  The  current  is  not  very  different  even  in  that  case.  We  also  looked  at  a  possible  effect  of 
secondary  electron  emission  (next  section). 
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10  Q 

*:  zero  seed  density  means  10  m'  in  the  near  vicinity  of  the  exposed  electrode  and  zero 
otherwise 

Figure  10a:  Current  on  the  lower  electrode  for  different  seed  electron  density 


Figure  10b:  Current  on  the  lower  electrode  for  different  seed  electron  density 


9 


Grant  FA8655-06-C-4004 

2nd  report  phase  3  -  January  2010 


3 

x  10'3 

. 

2.5 

^  " 

o 

o 

3 

1 

U) 

2 

10°  m'3 
-10-10m-3 

- 

1.5 

0  m'3  * 

- 

1 

0.5 

^  i  i  i  i  ,  1 _ i _ i _ i  .  i  i  .  ■  1 _ i _ i _ i _ ■  .ml _ i _ i _ l 

i  i  i  i  i  1 _ i _ i _ i _ 

io"9  io"8  io"7  io"6  io"5  10 


-4 


Time  (ns) 


Figure  11:  Total  energy  deposition  in  the  neutral  gas  for  different  seed  electron  density 


1.4  SENSITIVITY  TO  THE  ELECTRON  SECONDARY  EMISSION 

The  secondary  electron  emission  from  the  surface  due  to  ion  impact  is  another  key  parameter 
that  could  play  a  role  in  the  propagation  of  the  discharge.  However  as  shown  on  Fig.  12  the 
discharge  current  and  consequently  the  heating  properties  of  the  discharge  are  not  altered  by  the 
secondary  emission  coefficient. 
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Figure  12:  Current  on  the  lower  electrode  for  different  values  of  the  secondary  electron  emission 
coefficient  on  the  dielectric  barrier  surface 
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The  results  above  (only  slight  dependence  of  the  results  on  the  initial  seed  electron  density  or  on 
the  secondary  electron  emission)  are  very  intriguing  and  pose  the  question  of  the  possible  effects  of 
numerical  diffusion  (a  cathode  streamer  should  not  propagate  without  seed  electrons  or  when  the 
secondary  emission  coefficient  is  zero).  In  these  conditions  of  large  overvoltage  and  fast  rise  time 
numerical  diffusion  is  possible  in  spite  of  the  use  of  a  second  order  MUSCL  scheme  in  the  model 
and  this  point  needs  further  investigations. 

1.5  SENSITIVITY  TO  THE  AIR  FLOW  VELOCITY 

The  purpose  of  this  section  is  to  give  some  qualitative  preliminary  ideas  about  how  the  discharge 
and  its  effects  are  affected  by  an  external  airflow.  In  order  to  achieve  this,  new  boundary  conditions 
of  the  characteristic  type  have  been  implemented  for  subsonic  inflow  and  non  reflective  outflow  for 
the  Navier-Stokes  equations.  The  accuracy  of  the  external  airflow  computation  is  relatively  crude 
because  the  cell  number  has  to  be  limited  in  order  to  be  able  to  compute  the  discharge  in  a 
reasonable  time.  Furthermore  there  is  no  turbulence  model  added  to  the  Navier-Stokes  equations,  so 
we  expect  the  results  to  be  less  precise  as  the  velocity  increases.  However  the  aim  of  this  parametric 
study  is  only  to  draw  trends.  Figure  13  shows  the  external  airflow  obtained  for  50  m/s  incoming 
velocity.  Because  the  electrodes  are  located  at  3.5  cm  from  the  leading  edge  of  the  plate  and 
because  the  discharge  occurs  very  close  to  the  surface,  the  plasma  formation  is  not  so  much 
influenced  by  the  airflow.  Still  there  is  some  effect  as  shown  on  Fig.  14,  an  increase  in  the  external 
airflow  velocity  leads  to  a  lower  discharge  current  and  corresponding  heating  (see  Fig.  15)  and 
maximum  temperature  (see  Fig.  16).  Additional  investigation  are  required,  but  this  effect  is  believed 
to  be  higher  if  the  electrodes  are  placed  closer  to  the  leading  edge. 

Figure  17  shows  how  the  pressure  wave  generated  near  the  electrodes  interacts  with  the  external 
airflow.  Due  to  the  position  of  the  discharge  deep  into  the  boundary  layer,  the  interaction  is  low  at 
the  first  few  microseconds,  but  once  the  wave  has  gained  some  altitude,  it  encounters  the  stronger 
flow  and  begins  to  be  distorted. 
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Figure  13:external  airflow 
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4:  Current  on  the  lower  electrode  vs  time  for  different  velocities  of  the  external  flow 
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Figure  15:  maximum  air  temperature  reached 
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Figure  16:  Total  energy  deposition  in  the  neutral  gas 
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Figure  17:  Pressure  wave  interaction  with  the  external  airflow  showing  the  modification  of  the 
pressure  wave  (see  previous  report  and  Ref  [2])  at  two  different  times  (indicated  on  top  of  the 
figures)  ,  and  for  external  flows  from  0  to  150  m/s  (indicated  on  the  left  of  the  figures) 

1.6  CONCLUSION 

The  sensitivity  of  the  thermal  effects  created  by  a  nanosecond  pulsed  SDBD  on  the  neutral  gas 
has  been  investigated.  The  amplitude  of  the  voltage  pulse  seems  to  be  a  natural  and  efficient  driving 
parameter  of  the  effect.  The  duration  of  the  pulse  can  be  optimized  as  well.  The  rise  time  is  also  an 
important  parameter  that  must  be  investigated.  Some  trends  of  the  behavior  of  the  device  operating 
in  an  external  airflow  have  been  identified. 
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II.  ACCURACY  OF  THE  RESULTS:  A  MODEL  COMPARISON 

The  purpose  of  this  section  is  to  assess  the  sensitivity  of  the  results  to  the  numerical  schemes 
used  for  solving  the  charged  particles  transport  equations.  Two  numerical  methods  are  considered. 
The  first  one  is  based  on  the  exponential  scheme  of  Scharfetter  and  Gummel,  semi-implicitely 
coupled  with  Poisson’s  equation.  This  method  first  order  in  time  and  space.  The  second  method  is 
based  on  the  MUSCL  scheme  which  is  asynchronously  integrated  in  time  with  second  order  space 
accuracy  and  first  order  time  accuracy. 


L 


Figure  18:  Simulation  domain  for  model  comparison 

II.1  SINUSOIDAL  VOLTAGE  APPLIED  ON  THE  SDBD 

In  this  case  the  geometry  is  the  following:  L=2.8  mm,  w=0.1  mm,  er=4  and  the  upper  electrode 
length  is  0.35  mm.  The  mesh  is  composed  of  a  400x140  grid  of  uniform  square  cells  of  7  pm. 

A  8kV/8kHz  sinus  is  imposed  on  the  upper  electrode  for  two  periods.  Figure  19  &  20  shows  the 
current  response  of  each  simulation.  In  the  potential  rise  phase,  the  semi-implicit  scheme  solution 
consists  in  less  numerous  current  pulses  than  the  asynchronous  scheme  with  higher  peak  values. 
During  the  potential  drop  phase,  the  situation  is  inverted.  These  data  are  difficult  to  compare  with 
experiments  because  in  real  life  the  problem  is  not  homogenous  in  the  spanwise  direction.  However 
the  real  question  is  how  the  numerical  scheme  affects  the  computation  of  the  averaged  EHD  force. 
The  total  time-averaged  horizontal  EHD  force,  which  is  the  useful  parameter  if  ones  intends  to 
study  wall  jet  airflow  generation,  is  quite  of  the  same  order  of  magnitude  for  both  schemes.  The 
semi-implicit  scheme  gives  0.0212  N/m  whereas  the  asynchronous  scheme  gives  0.0236  N/m.  Some 
differences  are  present  in  the  space  distribution  of  the  force  as  shown  on  Figs.  21.  The  force  seems 
to  be  more  concentrated  in  the  electrode  tip  vicinity  for  the  asynchronous  scheme. 
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Figure  19:  Current  vs  voltage  for  the  semi-implicit  e  scheme 


Figure  20:  Current  vs  voltage  for  the  asynchronous  scheme 
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Figure  21:  Horizontal  time-averaged  EHD  force  for  semi-implicit  (top)  and  asynchronous  (bottom) 
schemes 


II.2  CONCLUSION 

Two  numerical  schemes  have  been  compared  in  terms  of  EHD  force  computation.  Although  the 
calculated  currents  are  not  identical  (number  of  pulses  per  cycle  and  values  of  the  peak  currents  are 
different),  the  trends  are  very  similar.  The  integrated  EHD  forces  are  equal  within  10%  and  the 
space  distributions  of  the  forces  are  similar. 
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II  SUMMARY  AND  CONLUSIONS  OF  THIS  WORK  ON  “STUDIES  OF  THE 
ELECTROHYDRODYNAMIC  FORCE  PRODUCED  IN  A  DIELECTRIC  BARRIER 
DISCHARGE  FOR  FLOW  CONTROL” 

III.l  ION  WIND 

The  mechanisms  of  ion  wind  generation  in  surface  dielectric  barrier  discharges  (SDBDs)  is 
identical  to  the  one  in  corona  discharges.  The  electrohydrodynamic  force  (EHD)  generated  in 
SDBDs  for  moderate  increase  rates  of  the  applied  voltage  is,  as  in  corona  discharges,  due  to  the 
development  of  a  unipolar  region  above  the  dielectric  surface,  where  ion  space  charge  is  dominant 
and  ions  transfer  momentum  to  the  neutral  gas.  For  very  fast  rise  (few  ns)  of  the  applied  voltage,  ion 
wind  is  no  longer  important  and  the  action  of  the  discharge  on  the  flow  is  different  (see  second  part 
of  this  conclusion).  SDBDs  are  more  interesting  for  applications  than  DC  corona  because  high 
current  breakdown  and  transition  to  the  arc  regime  are  prevented  by  the  use  of  a  dielectric  layer. 
There  is  no  evidence  that  SDBDs  can  generate  larger  ion  wind  than  DC  corona  discharges. 

Due  to  the  presence  of  a  dielectric  layer,  SDBDs  must  be  operated  with  a  time  varying  voltage. 
For  a  sinusoidal  regime  and  in  the  usual  electrode  configuration  of  SDBDs  for  flow  control 
applications,  transient  ion  clouds  develop  above  the  dielectric  surface  at  each  half  cycle.  During  the 
part  of  the  cycle  when  the  exposed  electrode  is  an  anode  (positive  part  of  the  cycle),  a  positive  ion 
cloud  forms  and  transfers  momentum  to  the  neutral  molecules.  A  negative  ion  space  charge  forms 
when  the  top  electrode  is  a  cathode  (negative  part  of  the  cycle)  and  the  EHD  force  is  due  to  negative 
ions  in  that  case.  The  contribution  of  positive  ions  or  negative  ions  to  the  overall  EHD  force 
depends  on  the  voltage  and  frequency  of  the  applied  voltage.  In  the  positive  regime,  the  formation 
of  ion  clouds  is  continuously  interrupted  by  high  current  pulses  (surface  streamers)  that  are  not 
efficient  for  ion  wind  production.  In  the  negative  regime  the  discharge  consists  of  high  frequency, 
low  amplitude  current  pulses  during  which  the  negative  ion  cloud  grows  continuously.  In  the 
negative  regime  a  large  force  per  unit  volume  also  exists  in  the  opposite  direction;  this  force  is 
however  limited  to  a  very  small  volume  corresponding  the  positive  ion  sheath  around  the  exposed 
electrode  (cathode  in  that  case) 

The  maximum  total  (space  and  time  integrated)  EHD  force  per  unit  length  of  the  electrode  that 
we  have  obtained  in  the  simulations  is  on  the  order  of  100  mN/m  (see  Fig.  22).  This  force  is 
distributed  along  the  surface  over  a  length  that  increases  with  the  amplitude  of  the  applied  voltage. 

The  best  discharge  efficiency  for  the  generation  of  the  EHD  force  is  obtained  in  the  simulations 
at  low  frequencies  and  high  voltages  as  seen  in  Fig.  22.  This  figure  shows  the  calculated  force  per 
unit  length  as  a  function  of  dissipated  power,  obtained  by  varying  the  amplitude  and  frequency  of 
the  applied  voltage.  Figure  22  shows  that  the  force  increases  as  a  function  of  power  for  all 
frequencies,  but  with  larger  slopes  at  low  frequencies.  At  “high  frequencies”  (eg  10  kHz),  the  force 
tends  to  saturate  at  lower  values  when  the  dissipated  power  increases.  The  relative  importance  of 
the  contribution  of  the  positive  and  negative  ions  is  also  shown  on  Fig.  22  (see  caption  to  Fig.  22). 
For  a  given  frequency,  the  relative  contribution  of  negative  ions  is  dominant  at  higher  voltages.  The 
transition  between  dominant  contribution  of  positive  and  negative  ions  occurs  at  lower  power  for 
higher  frequencies. 

The  model  results  show  excellent  qualitative  and  good  quantitative  agreement  with  experiments. 
Systematic  calculations  of  the  velocity  fields  corresponding  to  the  EHD  force  distributions  deduced 
from  this  plasma  model  would  be  useful.  The  model  results  can  also  be  used  to  provide  scaling  laws 
for  the  space  distribution  of  the  EHD  force  as  a  function  of  voltage  amplitude  and  frequency.  This 
would  provide  a  good  alternative  to  the  model  of  Susen  (eg  Susen  and  Huang,  44th  AIAA 
Aerospace  Sciences  Meeting  and  Exhibit,  9-12  January  2006,  Reno,  Nevada,  paper  AIAA  2006- 
877),  which  is  used  by  several  groups  in  the  US  but  is  lacking  of  solid  physical  basis. 
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Figure  22:  Calculated  EHD  force  parallel  to  the  surface,  per  unit  length  of  the  electrode,  as  a 
function  of  the  power  per  unit  length,  for  different  voltage  amplitudes  and  frequencies  of  the 
sinusoidal  voltage.  The  curves  are  plotted  for  constant  frequencies.  The  value  of  the  voltage 
amplitude  is  also  indicated  for  each  point.  The  X  on  the2,  5,  and  10  kHz  curves  shows  the  limit, 
along  the  curves,  where  the  contributions  of  positive  and  negative  ions  to  the  overall  force  are 
equal.  On  the  right  of  these  points  (larger  power)  the  contribution  of  negative  ions  is  dominant, 
while  the  contribution  of  positive  ions  is  dominant  on  the  left  of  the  points  marked  with  an  X.  The 
dielectric  layer  thickness  is  1  mm  and  the  length  and  height  of  the  simulation  domain  are  8  mm  and 
4  mm  respectively  (same  data  as  Fig.  14  of  previous  report). 

III.2  NANOSECOND  DISCHARGES 

For  a  given  power  dissipated  in  the  discharge,  Figure  22  shows  that  the  EHD  force  decreases 
when  the  frequency  increases.  This  trend  is  also  valid  for  other  voltage  waveforms.  When  using 
voltage  pulses,  the  EHD  force  decreases  when  the  voltage  rise  becomes  very  fast.  However  it  has 
been  shown  experimentally  by  Starikovskii  et  al.  that  high  voltage  pulses  with  very  short  rise  times 
(“nanosecond  voltage  pulses”  or  “nanosecond  regime”)  can  have  an  effect  on  the  flow. 

The  nanosecond  discharge  regime  has  been  simulated  in  this  work,  and  its  possible  aerodynamic 
effects  have  been  studied  by  coupling  the  discharge  model  with  Navier  Stokes  equations.  The 
results  confirm  that  the  EHD  force  generated  in  the  nanosecond  regime  is  negligible  and  show  that 
the  fast  heating  of  the  gas  in  the  vicinity  of  the  exposed  electrode  can  lead  to  a  large  temperature 
increase  in  a  short  time,  giving  rise  to  the  development  of  micro  shockwaves.  These  pressure  waves 
may  be  responsible  for  the  observed  aerodynamic  effects.  Further  work  is  needed  to  fully 
understand  the  interaction  of  the  pressure  waved  generated  by  the  nanosecond  discharges  with  the 
external  flow. 


18 


Grant  FA8655-06-C-4004 

2nd  report  phase  3  -  January  20 1 0 

REFERENCES 

[1]  A.  Starikovskii  ,  D.  Roupassov,  A.  Nikipelov,  and  M.  Nudnova  ,  SDBD  plasma  actuator  with 

nanosecond  pulse  periodic  discharge.  Plasma  Sources  Sci.  Technol.  18  034015  (2009) 

[2]  T.  Unfer,  and  J.P.  Boeuf  ,  Modeling  of  a  nanosecond  surface  discharge  actuator,  J.  Phys.  D: 
Appl.  Phys.  42  194017(2009) 


ATTACHED  PAPERS  PUBLISHED  IN  THE  FRAME  OF  THIS  CONTRACT 


Y  Lagmich,  Th  Callegari,  L  C  Pitchford  and  J  P  Boeuf,  Model  description  of  surface  dielectric 
barrier  discharges  for  flow  control ,  J.  Phys.  D:  Appl.  Phys.  41  095205  (2008) 


J.P.  Boeuf,  Y.  Lagmich,  L.  Pitchford,  Contribution  of  positive  and  negative  ions  to  the 
electrohydrodynamic  force  in  a  dielectric  barrier  discharge  plasma  actuator  operating  in  air, 
J.  Appl.  Phys.  106  023115  (2009) 


T.  Unfer,  and  J.P.  Boeuf ,  Modeling  of  a  nanosecond  surface  discharge  actuator,  J.  Phys.  D:  Appl. 
Phys.  42  194017(2009) 


19 


IOP  Publishing 


Journal  of  Physics  D:  Applied  Physics 


J.  Phys.  D:  Appl.  Phys.  41  (2008)  095205  (lOpp)  doi:l 0.1 088/0022-3727/4 1/9/095205 

Model  description  of  surface  dielectric 
barrier  discharges  for  flow  control 

Y  Lagmich,  Th  Callegari,  L  C  Pitchford  and  J  P  Boeuf 

LAPLACE,  Universite  de  Toulouse,  CNRS,  118  route  de  Narbonne,  31062  Toulouse,  France 

Received  29  January  2008,  in  final  form  29  February  2008 

Published  3  April  2008 

Online  at  stacks. iop.org/JPhysD/4 1/095205 

Abstract 

This  paper  presents  a  study  of  the  development  of  a  surface  dielectric  barrier  discharge  in  air 
under  conditions  similar  to  those  of  plasma  actuators  for  flow  control.  The  study  is  based  on 
results  from  a  2D  fluid  model  of  the  discharge  in  air  that  provides  the  space  and  time  evolution 
of  the  charged  particle  densities,  electric  field  and  surface  charges.  The  electrohydrodynamic 
(EHD)  force  associated  with  the  momentum  transfer  from  charged  particles  to  neutral 
molecules  in  the  volume  above  the  dielectric  layer  is  also  deduced  from  the  model.  Results 
show  that  the  EHD  force  is  important  not  only  during  the  positive  part  of  the  sinusoidal  voltage 
cycle  (i.e.  when  the  electrode  on  top  of  the  dielectric  layer  plays  the  role  of  the  anode)  but  also 
during  the  negative  part  of  the  cycle  (cathode  on  top  of  the  dielectric  layer).  During  the 
positive  part  of  the  cycle,  the  EHD  force  is  due  to  the  formation  of  a  positive  ion  cloud  that  is 
periodically  interrupted  by  high  current  breakdown.  The  EHD  force  during  the  negative  part  of 
the  cycle  is  due  to  the  development  of  a  negative  ion  cloud  that  continuously  grows  during  the 
successive  high  frequency  current  pulses  that  form  in  this  regime. 

(Some  figures  in  this  article  are  in  colour  only  in  the  electronic  version) 


1.  Introduction 

Surface  dielectric  barrier  discharges  (DBDs)  at  atmospheric 
pressure  can  generate  a  flow  or  modify  the  boundary  layer 
of  a  flow  and  have  been  proposed  as  actuators  for  flow 
control  [1-4].  The  momentum  transfer  from  charged  particles 
to  neutral  molecules  generates  an  electrohydrodynamic  (EHD) 
force  that  can  be  used  to  modify  the  airflow  profile  within 
the  boundary  layer  in  order  to  control  the  laminar-turbulent 
transition,  reduce  the  drag  and  reattach  or  stabilize  the  flow. 
The  ion  wind  in  corona  discharges  [5, 6]  is  a  good  example  of 
flow  generation  due  to  the  momentum  transfer  from  charged 
particles  to  neutral  molecules  in  a  gas  discharge. 

In  recent  papers  [7-9],  we  presented  studies  of  surface 
DBDs  in  pure  nitrogen  and  showed  that  the  EHD  force  was 
due  to  the  momentum  transfer  from  positive  ions  to  neutral 
molecules  during  the  formation  of  a  non-neutral  positive  ion 
cloud  above  the  surface  when  the  electrode  above  the  dielectric 
layer  plays  the  role  of  an  anode.  Once  the  size  and  density 
of  the  positive  ion  cloud  reach  critical  values,  a  high  current 
breakdown  occurs,  characterized  by  the  development  of  a  high 
density  filamentary  plasma  along  the  dielectric  surface.  The 
discharge  during  the  positive  part  of  the  cycle  is  composed 
of  successive  phases  of  ion  cloud  formation  and  high  current 


breakdown.  The  frequency  of  the  high  current  pulses  increases 
when  the  slope  of  the  voltage  waveform  increases.  The  spatial 
extension  of  the  positive  ion  cloud  (and  of  the  EHD  force) 
is  limited  by  these  breakdown  events  and  increases  when  the 
slope  of  the  applied  voltage  decreases  or  when  the  capacitance 
of  the  dielectric  layer  decreases. 

In  this  paper  we  study  a  surface  DBD  in  air.  We  show  that 
the  discharge  development  during  the  positive  part  of  the  cycle 
is  the  same  as  in  pure  nitrogen  but  that  the  negative  half  cycle 
is  different  because  of  the  negative  ion  generation  in  air.  The 
negative  part  of  the  cycle  is  also  composed  of  current  pulses, 
but  with  a  frequency  much  larger  than  in  the  positive  case 
and  of  much  lower  intensity.  A  negative  ion  cloud  forms  and 
continuously  grows  during  the  negative  part  of  the  cycle.  This 
ion  cloud  can  generate  an  EHD  force  as  large  or  even  larger 
than  the  force  generated  by  positive  ions  during  the  positive 
phase. 

In  section  2  we  recall  the  principles  of  the  model 
and  of  the  force  generated  by  momentum  transfer  between 
charged  particles  and  neutrals.  Section  3  describes  the 
discharge  mechanisms  for  positive  and  negative  ramp  voltage 
waveforms.  Conclusions  are  drawn  in  section  4. 
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2.  EHD  force  per  unit  volume  and  discharge  model 

As  discussed  in  previous  publications  [7-9],  the  EHD  force 
per  unit  volume  in  electric  discharges  is  due  to  momentum 
transfer  from  charged  particles  to  neutral  particles.  Neglecting 
the  mean  velocity  of  the  neutral  particles  with  respect  to  the 
charged  particle  mean  velocities,  the  EHD  force  can  be  written 
as  (the  indices  ‘e’,  T  and  ‘n’  referring  to  electrons,  positive 
ions  and  negative  ions,  respectively) 

/e,  i,n  =  ^e,i,n^e,i,n^me,i,n^e,i,n?  (1) 

where  are  the  charged  particle  densities,  me,i,n  their 
masses,  vme,i,n  the  momentum  exchange  frequencies  for 
electron-neutral,  positive  ion-neutral  and  negative  ion-neutral 
collisions  and  ue^n  the  charged  particle  mean  velocities. 

Writing  the  mobility  of  species  V,  as  fis  =  e/(msvms), 
the  total  force  per  unit  volume/ is 


j  —  •  v 

Ml  Me  Mn 

In  the  conditions  of  atmospheric  discharge  plasma  actuators, 
we  have  seen  (see  [7-9])  that  the  dominant  contribution  to  the 
EHD  force  comes  from  the  drift  terms  of  the  current  densities 
of  equation  (2),  i.e. 

js  «  ensfisE.  (3) 

We  can  therefore  write 

/=  -  -  —  -  —  *  e(n,  -  ne  -  na)E.  (4) 

Mi  Me  Mn 

The  force  per  unit  volume  acting  on  the  neutral  molecules  is 
therefore  equal  to  the  Coulomb  force  acting  on  the  charged 
particles,  which  means  that  the  momentum  gained  by  the 
charged  particles  in  the  electric  field  is  exactly  and  locally 
balanced  by  collisions,  and  entirely  transmitted  to  neutral 
molecules.  One  consequence  of  equation  (4)  is  that  the  EHD 
force  is  zero  in  a  first  approximation  for  a  quasi-neutral  plasma 
(see  [7]  for  a  more  explicit  discussion).  Equation  (4)  indicates 
that  the  EHD  force  is  large  in  regions  of  the  discharge  with 
a  large  electric  field  and  where  the  space  charge  is  non-zero 
(positive  or  negative).  This  is  the  case,  for  example,  in  the 
cathode  sheath  regions  of  glow  discharges  or  in  the  drift  region 
of  corona  discharges. 

In  this  paper,  the  space  and  time  variations  of  the  EHD 
force  are  deduced  from  a  numerical  model  of  the  discharge. 
The  model  is  based  on  fluid  equations  for  electrons  and  ions, 
coupled  with  Poisson’s  equation  for  the  electric  field.  The 
time  dependent  electron  and  ion  continuity  equations  with  a 
drift-diffusion  flux  are  coupled  with  Poissons ’s  equation  and 
integrated  in  time  using  the  Scharfetter-Gummel  scheme  with 
a  semi-implicit  method  (see  [7-9]  and  references  therein). 

In  this  study,  a  very  simple  model  of  the  plasma  chemistry 
has  been  considered,  with  one  type  of  positive  ion  and  one 
type  of  negative  ion  species  (for  air)  and  no  complex  plasma 
chemistry  is  included.  Negative  ions  in  air  are  supposed  to  be 
formed  in  two-body  and  three-body  reactions.  Electron-ion 


and  ion-ion  recombination  are  taken  into  account  but  negative 
ion  detachment  is  not  considered. 

The  continuity  equation  for  the  charged  particles  of  type 
s  (s  =  e,  i,  or  n)  reads 

dns 

-^  +  V.Ts  =  Ss,  (5) 

ot 

where  the  charged  particle  flux  Ts  has  a  drift-diffusion  form 

Ts  =  ±nsfisE  -  DsWns ,  (6) 

where  the  +  sign  applies  to  positive  ions  and  the  —  sign  is  for 
electrons  and  negative  ions.  fis  and  Ds  are  the  mobility  and 
the  diffusion  coefficient  of  species  s . 

The  source  terms  Ss  of  the  continuity  equations  are 
defined  by 

SG  =  (a-  ?7)|re|  - 
=  or|re|  —  reinQrii  -  rmnxn^ 

Sn  =  rj\VQ\-  rinmnn ,  (7) 

a  and  r)  are  the  ionization  and  attachment  coefficient, 
respectively;  re i  and  rm  are  the  electron-positive  ion  and 
positive  ion-negative  ion  recombination  coefficients. 

The  continuity  equations  (5)  with  the  flux  from  the 
momentum  equations  (6)  and  the  source  terms  from 
equation  (7)  are  coupled  with  Poisson’s  equation: 

V  -  (eE)  =  —(rii  -ne-  nn)  +  —  cr8s  (8) 

£o  s  o 

a,  in  the  last  term  of  equation  (8),  is  the  surface  charge  density, 
which  is  non-zero  only  on  the  dielectric  layer  surface,  as 
expressed  by  the  Dirac  function  8S.  The  surface  charge  density 
is  calculated  self-consistently  by  integrating  the  electron  and 
ion  fluxes  to  the  surface.  Electrons  and  ions  in  each  surface 
element  are  supposed  to  recombine  instantly  with  oppositely 
charged  particles  on  the  surface  if  present. 

The  current  calculated  in  the  model  and  displayed  in  some 
of  the  figures  below  is  the  displacement  current  on  the  bottom 
electrode.  Because  of  the  boundary  conditions  (zero  current 
out  of  the  simulation  domain),  this  current  is  also  equal  to 
the  sum  of  the  electron,  ion  and  displacement  currents  on 
the  electrode  above  the  surface  (the  displacement  current  is 
calculated  on  both  sides  of  the  electrode). 

When  the  electric  field  is  directed  towards  the  dielectric 
layer  surface,  secondary  electron  emission  is  taken  into  account 
by  the  boundary  condition:  7e,±  =  yJi,±  where  /e,±  and  J\  j_ 
are  the  components  of  the  electron  and  ion  current  densities 
perpendicular  to  the  dielectric  surface  and  y  is  the  secondary 
electron  emission  coefficient  due  to  ion  impact,  y  is  set  to  0.05 
in  the  calculations  below. 

The  simulation  domain  is  shown  in  figure  1 .  The  dielectric 
layer  permittivity  sr  is  set  to  5,  unless  otherwise  indicated.  The 
applied  voltage  waveform  between  the  electrodes  is  supposed 
to  be  linearly  increasing  with  time  (ramp  voltage),  with  a 
slope  rjv,  on  the  order  of  100  V  /xs-1  or  so  (this  is  consistent 
with  typical  experiments  for  sinusoidal  voltage  waveform  of 
5-20  kV  amplitude  and  frequencies  in  the  1-10  kHz  range). 
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Figure  1.  Simulation  domain.  The  calculations  below  have  been 
performed  mainly  for  the  two  sets  of  conditions:  (L  =  4  mm, 
w  =  0.5  mm,  h  —  1.5  mm,  /  =  0.25  mm)  and  (L  =  8  mm, 
w  =  1  mm,  h  —  3  mm  /  =  0.5  mm). 
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Figure  2.  ( a )  Reduced  ionization  and  attachment  coefficients 
used  in  the  model  as  a  function  of  reduced  electric  field  (at  300  K); 

0 b )  reduced  electron  and  ion  mobilities  as  a  function  of  reduced 
electric  field  (at  300  K).  The  air  data  are  used  in  the  present  model; 
the  nitrogen  data  also  shown  for  comparison  in  this  figure  were  used 
in  [7-9]. 

The  ionization  and  attachment  coefficients  are  supposed 
to  depend  on  the  local  reduced  electric  field  as  in  [7-9]  and 
are  obtained  with  the  BOLSIG+  electron  Boltzmann  equation 
solver  [10]  (three-body  attachment  in  oxygen  is  properly  taken 
into  account  with  a  separate  cross  section).  The  electron¬ 
positive  ion  and  positive  ion-negative  ion  recombination 
coefficients  re i  and  rm  are  supposed  to  be  equal  to  2  x 
10_7cm3s_1.  The  electron  and  ion  diffusion  coefficients 
are  chosen  so  that  De//xe  =  IV  and  A,n//^i,n  =  0.01  V. 
Ionization  and  attachment  coefficients  and  mobilities  used  for 
nitrogen  in  [7-9]  and  in  air  in  this  study  are  shown  in  figure  2. 
A  constant  mobility  has  been  used  for  positive  and  negative 
ions  in  air. 

A  similar  physical  model  of  surface  DBD  plasma  actuators 
has  been  recently  published  by  Likhanskii  et  al  [1 1]. 


Figure  3.  Calculated  current  waveform  in  air  in  the  geometry  of 
figure  1  and  for  a  positive  voltage  slope  (also  shown) 

T]V  =  +300  V  /xs_1.  The  symbols  indicate  the  times  at  which  the 
charged  particle  densities  and  electric  potential  are  shown  in 
figure  4.  The  permittivity  of  the  dielectric  layer  er  is  equal  to  5.  The 
dimensions  of  the  simulation  domain  (see  figure  1)  are  (L  =  4  mm, 
w  =  0.5  mm,  h—  1.5  mm,  /  =  0.25  mm).  The  simulation  domain  is 
200  x  400,  i.e.  the  grid  spacing  is  10  /xm. 

Finally  note  that  the  numerical  method  used  in  this  paper  to 
solve  equations  (5)-(8)  is  first  order  accurate.  More  accurate, 
2nd  order  explicit  methods  are  being  developed  [12]  to  validate 
the  results  obtained  with  this  model. 

3.  Surface  DBD  in  air  in  the  plasma  actuator 
geometry 

The  results  presented  here  have  been  obtained  in  air  with  a 
simplified  model  taking  into  account  only  one  type  of  positive 
ion  and  one  type  of  negative  ion  (see  above). 

Results  are  presented  for  a  positive  ramp  voltage  (‘positive 
discharge’),  i.e.  linearly  increasing  positive  voltage  applied  to 
the  electrode  above  the  dielectric  layer,  the  electrode  below 
the  dielectric  layer  being  grounded  (the  electrode  below  the 
dielectric  layer  is  the  cathode),  and  for  a  negative  ramp  voltage 
(‘negative  discharge’),  i.e.  linearly  increasing  negative  voltage 
applied  to  the  electrode  above  the  dielectric  layer,  the  electrode 
below  the  dielectric  layer  being  grounded  (the  electrode  above 
the  dielectric  layer  is  the  cathode). 

We  find  that,  for  positive  ramp  voltage  waveform  (anode 
above  the  dielectric  surface),  most  of  the  features  described  in 
previous  papers  [7-9]  for  nitrogen  are  qualitatively  valid  for 
air.  For  a  negative  ramp  voltage  waveform  (cathode  above  the 
dielectric  surface)  a  significant  EHD  force,  directed  away  from 
the  top  electrode  (as  in  the  positive  case),  is  present  above  the 
surface  and  is  due  to  the  formation  of  a  negative  ion  cloud. 
Note  that  in  pure  nitrogen,  we  found  in  [8]  that  the  EHD  force 
for  a  negative  ramp  voltage  was  negligibly  small. 

We  present  results  in  air  for  positive  and  negative  ramp 
voltage  waveforms  in  the  two  subsections  below. 

3.1.  Positive  ramp  voltage  waveform 

3.1.1.  Time  evolution  of  the  positive  surface  DBD.  The 
calculated  current  waveform  in  air  shown  in  figure  3  for  a 
positive  ramp  voltage  (anode  above  the  dielectric  surface) 
consists  of  high  intensity  pulses  separated  by  low  current 
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Figure  4.  Distribution  of  charged  particle  densities  and  potential  at  different  times  for  a  discharge  in  air  in  the  conditions  of  figure  3  for  a 
ramp  voltage  r/y  =  +300  Y  /x s-1.  The  colour/greytone  scale  for  the  charged  particle  densities  is  a  log  scale  over  four  decades.  The  voltage 
contours  are  indicated  in  kilovolts.  (Colour  online  only.) 


phases,  and  is  similar  to  the  current  waveform  calculated  in 
nitrogen  (see  [7-9]). 

We  see  in  figure  4  that  the  low  current  phase  is  associated 
with  the  formation  of  a  positive  ion  cloud  above  the  dielectric 
surface.  This  ion  cloud  forms  and  starts  to  grow  as  soon 
as  the  potential  drop  along  the  dielectric  surface  becomes 
sufficient  for  the  discharge  to  become  self-sustained.  Under 
the  assumptions  of  our  model,  electrons  are  generated  at  the 
dielectric  surface  by  ion  bombardment.  They  multiply  in  the 
large  electric  field  above  the  surface.  In  the  low  current  phase 
(e.g.  times  t\  and  ^  in  figure  4),  the  field  distortion  due  to  the 
ion  space  charge  build  up  is  not  sufficient  to  lead  to  a  significant 
growth  of  the  electron  density,  and  the  electron  density  is 
small  with  respect  to  the  positive  ion  density  (no  plasma). 
Due  to  the  continuous  increase  in  the  applied  voltage  and  to 
the  dielectric  nature  of  the  surface,  the  ion  cloud  grows  and 
expands  along  the  surface.  When  the  ion  space  charge  reaches 
a  critical  value,  around  time  t$  =  25.97  /xs,  the  geometric  field 
becomes  significantly  distorted.  The  reconfiguration  of  the 
electric  field  due  to  the  plasma  formation  leads  to  a  significant 
increase  in  the  electron  multiplication  and  of  the  discharge 
current  (current  peaks  in  figure  3).  A  quasi-neutral  plasma  thus 
forms  and  a  filament  propagates  with  a  high  velocity  along  the 
surface  (times  1 3,  £4  of  figure  4).  At  the  end  of  the  current 
pulse  (around  time  *5)  most  of  the  electrons  are  attached  and 
the  post-discharge  plasma  is  composed  of  mainly  positive  and 
negative  ions.  The  dielectric  layer  surface  is  positively  charged 
and  the  electric  potential  immediately  above  the  surface  is 
near  the  anode  potential.  After  time  t$,  the  post-discharge 
plasma  diffuses  and  recombines.  Since  the  applied  voltage  is 
supposed  to  increase  linearly,  the  same  phenomena  (formation 
and  growth  of  a  positive  ion  cloud,  followed  by  high  current 
breakdown)  are  repeated  at  a  later  time  (around  t  =  52  /xs,  see 
figure  3)  and  a  periodic  regime  is  reached  in  the  simulation. 

It  is  interesting  to  look  more  closely  at  the  time  evolution  of 
the  electric  field  above  the  dielectric  surface.  Figures  5(a)  and 


Figure  5.  Distribution  along  the  surface  (x  direction,  see  figure  1) 
of  the  electric  field  above  the  surface  at  different  times  (a)  during 
the  ion  cloud  expansion  (low  current  phase)  and  (b)  during  the  high 
current  breakdown.  Same  conditions  as  figures  3  and  4. 

(b)  show  the  distribution  of  the  total  electric  field  immediately 
above  the  surface,  as  a  function  of  position  along  the  surface, 
before  and  after  high  current  breakdown,  respectively.  At  time 
t  =  1.66 /xs  (figure  5(a))  the  electric  field  is  not  perturbed 
by  the  space  charge  and  is  only  defined  by  the  electrode  and 
dielectric  arrangement.  The  electric  field  distortion  due  to 


4 


J.  Phys.  D:  Appl.  Phys.  41  (2008)  095205 


Y  Lagmich  et  al 


the  ion  cloud  appears  quickly  (see  the  field  at  t  =  5.2  /xs  in 
figure  5(a)).  The  field  decreases  in  the  vicinity  of  the  electrode 
tip,  and  increases  in  the  ion  cloud,  away  from  the  electrode  tip. 
The  region  of  large  electric  field  extends  as  the  ion  cloud  grows, 
until  breakdown  occurs  around  time  t  =  25.96-25.97  /xs  (see 
figures  5(a)  and  (b)). 

The  field  evolution  above  the  surface  after  t  =  25.97  /xs  in 
figure  5(b)  is  very  similar  to  the  electric  field  during  streamer 
propagation.  A  plasma  channel,  characterized  by  a  drop  of  the 
electric  field,  forms.  This  channel  is  surrounded  by  two  large 
electric  field  regions  on  the  cathode  side  and  on  the  anode 
side,  which  can  be  qualified  as  cathode  and  anode  streamers. 
The  cathode  streamers  propagates  at  a  high  velocity  (about 
1.5  mm  in  7  ns  between  t  =  25.973  /xs  and  t  =  25.98  /xs,  i.e. 
2  x  107  cms-1).  Note  that  photoionization  is  not  included  in 
the  model  and  that  the  cathode  streamer  propagation  in  the 
simulations  is  due  to  secondary  emission  from  the  dielectric 
surface  and/or  to  charges  remaining  in  the  volume  over  the 
dielectric  surface  from  previous  current  pulses.  Since  the 
discharge  takes  place  very  close  to  the  surface  it  seems 
reasonable  to  assume  that  photoionization  does  not  play  an 
important  role  in  these  surface  DBDs  (even  a  very  low  electron 
emission  from  the  surface  can  provide  a  fast  propagation  of 
the  discharge  along  the  surface).  Nevertheless,  the  exact 
propagation  velocity  of  the  streamer  along  the  surface  and  the 
time  interval  between  successive  pulses  are  very  difficult  to 
predict  accurately.  Moreover  the  numerical  method  being  first 
order  in  space,  this  can  have  some  influence  of  the  propagation 
velocity  and  on  the  time  interval  between  current  pulses.  Work 
is  in  progress  to  perform  sensitivity  analyses  of  the  results  to 
parameters  such  as  the  secondary  emission  coefficient,  and  to 
study  the  effect  of  the  accuracy  of  the  numerical  method  on  the 
results,  by  using  second  order  methods  as  described  in  [12]. 
Preliminary  results  show  that  the  main  conclusions  of  this 
paper  are  not  altered  when  the  secondary  emission  coefficient 
is  varied  or  when  a  more  accurate  method  [12]  is  used.  These 
results  will  be  discussed  in  forthcoming  publications. 

3.1.2.  EHD  force.  As  we  found  in  the  case  of  a  positive 
ramp  voltage  in  pure  nitrogen  in  [8, 9],  the  EHD  force  in  air 
is  important  only  during  the  low  current  phases  between  high 
current  pulses  (see  figure  5).  Figure  6  shows  the  variations  as 
a  function  of  time  t ,  of  the  space  integrated  EHD  force  parallel 
to  the  surface,  averaged  in  the  time  interval  [0,t]  and  defined  as 

F\i(t)  =  ~  [*  Fn(t')dt'  =  -  fdt'  f/n(x,y,Odxdy, 

t  Jo  t  Jo  J 

(9) 

where  the  instantaneous  EHD  force  per  unit  volume,  /,  is 
defined  in  equation  (4).  F||  is  a  force  per  unit  length  of  the 
electrode. 

The  EHD  force  F||  in  figure  6  is  represented  for  different 
values  of  the  ramp  voltage  slope  and  for  different  combinations 
of  dielectric  layer  thickness  and  permittivity.  The  oscillations 
of  the  force  in  figure  6  are  associated  with  the  current  pulses 
(the  force  increases  only  during  the  low  current  phase  between 
the  pulses).  We  see  that  the  frequency  of  these  oscillations 
increases  with  increasing  voltage  slope,  and  with  the  dielectric 


Time  (ps) 

Figure  6.  Variations  as  a  function  of  time  t,  of  the  EHD  force  per 
unit  length  (averaged  in  time  between  0  and  t ,  see  equation  (9)) 
parallel  to  the  surface,  for  different  slopes  of  the  positive  ramp 
voltage  and  for  different  values  of  the  dielectric  thickness  and 
permittivity.  Air,  simulation  domain  of  dimensions  L  =  8  mm  and 
w  +  h  =  4  mm,  /  =  0.5  mm  (see  figure  1). 


Figure  7.  Extension  of  the  positive  ion  cloud  along  the  surface  at 
breakdown  for  a  positive  voltage  slope,  as  a  function  of  the 
parameter  r]Vsr/w,  for  different  values  of  the  dielectric  permittivity 
in  air  for  simulation  domain  dimensions  L  =  8  mm,  w  =  1  mm  and 
h  =  3  mm,  /  =  0.5  mm. 

permittivity  to  thickness  ratio  (related  to  the  capacitance  of  the 
dielectric  layer),  as  was  already  found  in  pure  nitrogen  [8]. 
However,  we  note  that  the  time  averaged  EHD  force  tends  to 
reach,  at  steady  state,  a  value  on  the  order  of  0.12mNm_1, 
which  is  not  very  sensitive  to  these  parameters.  We  also  find, 
as  in  [8,  9]  for  pure  nitrogen,  that  although  the  space  and 
time  integrated  force  is  not  very  dependent  on  the  conditions, 
the  space  distribution  of  the  force  over  the  surface  is  very 
sensitive  to  the  ramp  voltage  slope  and  to  the  dielectric 
capacitance. 

The  length  of  the  ion  cloud  just  before  breakdown,  i.e.  the 
spatial  extension  of  the  EHD  force,  is  represented  in  figure  7 
as  a  function  of  r]v£r/w,  which  is  proportional  to  the  charging 
current  density  of  the  dielectric  layer.  We  see  in  figure  7  that  the 
spatial  extension  of  the  EHD  force  along  the  surface  depends 
only  on  the  scaling  parameter  r]y£r/w,  as  was  already  shown 
in  the  case  of  pure  nitrogen  in  [9]. 

When  the  slope  of  the  ramp  voltage  increases,  or  when  the 
capacitance  of  the  dielectric  layer  increases,  the  size  of  the  ion 
cloud  at  breakdown,  and  thus  the  spatial  extension  of  the  EHD 
force,  decreases. 
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Figure  8.  Frequency  of  the  high  current  pulses  and  expansion 
velocity  of  the  ion  cloud  as  a  function  of  the  positive  voltage  slope 
for  two  values  of  the  permittivity  of  the  dielectric  layer  and  for 
simulation  domain  dimensions  L  =  8  mm,  w  =  1  mm  and 
h  —  3  mm,  l  =  0.5  mm. 


It  is  interesting  to  look  at  the  variations  of  the  expansion 
velocity  of  the  ion  cloud  as  a  function  of  the  voltage  slope 
rjv  shown  in  figure  8.  We  see  that  this  velocity  increases 
practically  linearly  with  r]V  and  is  on  the  order  of  a  few 
100  ms-1  to  lkms-1  under  typical  conditions.  The  scaling 
with  r\  v  sr  /  uj  is  not  as  clear  here,  and  we  have  plotted  separately 
in  figure  8  the  results  as  a  function  of  r]V  for  two  different 
values  of  the  permittivity.  The  ion  cloud  expansion  velocity 
increases  with  rjv  but  with  slightly  different  slopes  for  different 
values  of  the  permittivity  sr  (for  a  given  dielectric  layer  width). 
The  same  remark  is  valid  for  the  current  pulse  frequency 
which  is  also  represented  in  figure  8.  Since  the  expansion 
length  of  the  positive  ion  cloud  is  limited  by  breakdown, 
this  length  (shown  in  figure  7)  is  directly  proportional  to 
the  ratio  of  the  ion  velocity  to  the  pulse  frequency  shown  in 
figure  8. 
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Figure  9.  ( a )  Calculated  current  waveform  in  air  in  the  geometry 
of  figure  1  and  for  a  negative  voltage  slope  (also  shown) 
rjv  =  —300  V  /x s-1.  The  symbols  indicate  the  times  at  which  the 
charged  particle  densities  and  electric  potential  are  shown  in 
figure  11.  ( b )  Detailed  view  of  a  single  current  pulse.  The  symbols 
indicate  the  times  at  which  the  charged  particle  densities  and  electric 
potential  are  shown  in  figure  10.  The  dimensions  of  the  simulation 
domain  (see  figure  1)  are  L  =  4  mm,  w  =  0.5  mm,  h  =  1.5  mm  and 
/  =  0.25  mm.  The  permittivity  of  the  dielectric  layer  sr  is  equal  to  5. 
The  simulation  domain  is  200  x  400,  i.e.  the  grid  spacing  is  10  /xm. 


3.2.  Negative  ramp  voltage  waveform 

3.2.1.  Time  evolution  of  the  negative  surface  DBD.  For 
negative  ramp  voltages  (cathode  above  the  dielectric  layer),  the 
current  exhibits  higher  frequency  pulses  of  smaller  amplitudes, 
as  shown  in  figure  9(a).  The  frequency  of  the  oscillations  is 
about  1  MHz  in  the  conditions  of  figure  9,  and  the  duration  of 
the  current  pulse  is  on  the  order  of  100  ns  (see  figure  9(b)).  We 
found  similar  current  oscillations  in  pure  nitrogen  (see  [8]). 
Note  that  in  air  (figure  9)  the  minimum  current  in  the 
oscillations  increases  slightly  with  time. 

Figure  10  shows  the  space  distribution  of  the  electron 
density,  positive  and  negation  ion  densities  and  electric 
potential  at  five  different  times  of  a  single  pulse,  indicated 
in  figure  9(b).  We  see  that  on  the  time  scale  of  one  current 
pulse  the  electron  density  is  completely  modulated  while  the 
ion  densities  slightly  change  in  time  but  do  not  decay  to  zero. 
At  time  4  before  the  current  increase,  the  electron  density  is 
practically  zero.  We  see  then,  at  times  4  to  4,  an  increase 
in  the  electron  density  leading  to  the  formation  of  a  quasi¬ 
neutral  plasma  in  the  vicinity  of  the  electrode  tip.  The  electron 
density  extends  beyond  the  plasma  region  (times  4,  4),  in  a 
region  which  is  dominated  by  a  negative  space  charge.  The 


electron  density  then  decays  and  goes  to  zero  (time  4)  due 
to  the  charging  of  the  dielectric  surface  by  the  electrons, 
leading  to  a  potential  drop  along  the  dielectric  surface  and 
to  a  decrease  in  the  electron  multiplication  below  that  of  the 
self-sustaining  condition.  The  decay  of  the  electron  density  is 
also  due  to  electron  attachment,  but  the  current  oscillations  are 
not  due  to  attachment  as  it  is  in  the  Trichel  regime  of  negative 
corona  discharges.  Note  that  we  have  shown  that  the  current 
oscillations  in  the  negative  regime  also  exist  in  pure  nitrogen 
where  attachment  is  absent  [8].  We  note  that  at  the  end  of 
the  current  pulse  (time  4)  the  densities  of  positive  ions  and 
negative  ions  at  the  extremity  of  the  ion  cloud  away  from  the 
top  electrode  have  increased,  the  negative  ion  density  being 
however  larger  than  the  positive  ion  density.  This  is  because  of 
electron  impact  ionization  and  attachment  during  the  electron 
transport  above  the  surface  (times  4-4).  We  can  therefore 
conclude  that  the  ion  space  charge  above  the  dielectric  surface 
is  negative  on  the  average  and  that  it  expands  at  each  current 
pulse  due  to  electron  attachment  (and  negative  ion  drift).  We 
will  see  below  that  the  EHD  force  during  the  negative  voltage 
ramp  is  due  to  the  negative  ion  space  charge  (negative  ion 
cloud)  that  develops  above  the  surface  and  expands  at  each 
current  pulse. 
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Figure  10.  Distribution  of  charged  particle  densities  and  potential  at  different  times  for  a  discharge  in  air  in  the  conditions  of  figure  9  for  a 
ramp  voltage  r]V  =  —300  V  /xs_1,  over  one  single  current  pulse  (times  are  indicated  on  the  current  pulse  of  figure  9(b)).  The  colour/greytone 
scale  for  the  charged  particle  densities  is  a  log  scale  over  four  decades.  The  voltage  contours  are  indicated  in  kilovolts.  (Colour  online  only.) 
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Figure  11.  Distribution  of  the  density  of  positive  ions  nx,  of  negative  ions  nn ,  of  the  net  negative  charged  particle  density  n  =  nQ  +  nn  —  nx, 
and  of  the  electric  potential  at  different  times  for  a  discharge  in  air  in  the  conditions  of  figure  9,  over  several  current  pulses  (times  are 
indicated  on  the  current  plot  in  figure  9(a)).  Only  the  positive  values  of  n  are  represented  (n  is  negative  in  the  positive  ion  sheath  near  the 
electrode  tip).  Note  that  the  colour/greytone  scale  for  the  charged  particle  densities  is  a  log  scale  over  three  decades,  and  not  on  four 
decades,  as  in  figure  10.  The  voltage  contours  are  indicated  in  kilovolts.  (Colour  online  only.) 


Figure  11  shows  the  space  distributions  of  the  electron 
density,  positive  ion  density,  net  density  of  charges  and  elec¬ 
tric  potential  at  different  times,  but  on  a  longer  time  scale  than 
in  figure  10.  The  distributions  are  shown  at  the  times  indicated 
in  figure  9(a) ,  i.e.  at  the  maxima  of  every  fourth  current  pulse, 
starting  at  t\  —  5.88  /xs.  We  see  clearly  in  this  figure  that  the 
positive  and  negative  ion  clouds  expand  continuously  during 
the  successive  pulses,  and  that,  away  from  the  electrode  tip, 
the  space  charge  is  essentially  negative.  Due  to  the  charging 


of  the  dielectric  surface,  the  length  of  the  electron  paths  above 
the  electrode  increases  at  each  current  pulse,  and  this  is  associ¬ 
ated  with  the  expansion  of  the  positive  and  negative  ion  clouds 
above  the  surface.  Note  that  the  region  in  the  vicinity  of  the 
electrode  tip  (which  plays  the  role  of  the  cathode)  is  similar  to 
the  cathode  region  of  a  glow  discharge  (with  a  strong  field  dis¬ 
tortion  due  to  the  electrode  geometry,  as  in  a  corona  discharge). 

It  is  interesting  to  look  more  closely  at  the  time  evolution 
of  the  surface  charge  and  electric  field  above  the  surface. 
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Figure  12.  Space  distribution  of  the  surface  charge  in  the  x 
direction  (see  figure  1)  along  the  surface  at  different  times  during 
the  discharge  of  figures  9-11. 

Figure  12  shows  the  distribution  of  the  charge  along  the  surface 
(due  to  the  charging  by  electrons  and  negative  ions)  at  different 
times  considered  in  figure  1 1  (and  indicated  in  figure  9(a)).  We 
see  that  the  extension  of  the  negative  charges  along  the  surface 
increases  with  time,  pulse  after  pulse.  This  is  consistent  with 
the  results  of  figure  11  above  where  it  was  shown  that  the 
electron  density  above  the  dielectric  layer  extends  a  little  more 
at  each  current  pulse,  due  to  the  charging  of  the  surface.  The 
maximum  of  the  surface  charge  is  located  a  few  100  gm  away 
from  the  electrode  tip  and  does  not  move  significantly  during 
the  time  interval  considered  (about  20  /as).  This  maximum  of 
the  surface  charge  and  the  charge  along  the  surface  within 
the  ion  clouds  increase  practically  linearly  with  time,  and 
at  a  constant  rate  along  the  surface.  The  gradient  of  the 
surface  charge  after  the  maximum  is  almost  constant  in  space 
and  time. 

Figure  13  shows  the  distribution  of  the  electric  field 
immediately  above  the  surface,  at  the  same  times  as  in 
figures  1 1  and  12.  The  electric  field  is  large  near  the  electrode 
tip,  and  decreases  sharply  away  from  it.  We  see  that  the  electric 
field  presents  a  relative  maximum  that  moves  away  from  the 
electrode  tip.  This  maximum  electric  field  is  on  the  order 
of  the  breakdown  field  in  air  i.e.  30  kV  cm-1  (the  value  for 
which  ionization  balances  attachment).  The  evolution  of  the 
electric  field  above  the  surface  is  associated  with  the  charging 
of  the  dielectric  and  the  evolution  of  the  space  charge  above 
the  surface. 

3.2.2.  EHD  force.  As  mentioned  above,  the  EHD  force 
for  a  negative  voltage  in  air  can  be  large,  in  contrast  to 
the  case  of  pure  nitrogen  (see  [8]).  The  EHD  force  is 
due  to  the  development  of  a  negative  space  charge  above 
the  surface.  Note  that  the  expansion  of  the  negative  space 
above  the  surface  is  continuous  and  is  not  interrupted  by 
high  current  breakdown,  as  in  the  case  of  the  positive  voltage 
(see  section  3.1).  Therefore,  the  length  of  the  negative 
ion  space  charge  above  the  dielectric  surface  is  not  limited 
by  breakdown  but  rather  by  the  duration  of  the  negative 
phase,  i.e.  by  the  frequency  of  the  applied  voltage.  One  can 
therefore  expect  that  the  EHD  force  during  the  negative  part 


Figure  13.  Space  distribution  of  electric  field  above  the  surface  at 
different  times  during  the  discharge  of  figures  9-1 1 . 


Figure  14.  Time  variations  of  the  EHD  force  per  unit  length 
parallel  to  the  surface  for  different  slopes  of  the  negative  ramp 
voltage  and  for  different  values  of  the  dielectric  thickness  and 
permittivity,  in  air  and  for  simulation  domain  dimensions  L  =  8  mm 
and  w  +  h  =  4  mm,  /  =  0.5  mm  (see  figure  1). 

of  the  cycle  will  play  a  more  important  role  at  lower  voltage 
frequencies. 

Figure  14  shows  the  time  evolution  of  the  space  and  time 
averaged  EHD  force  per  unit  length  parallel  to  the  surface, 
F ||  (see  equation  (9)),  for  a  negative  ramp  voltage  in  different 
conditions  (different  dielectric  layer  permittivity  and  thickness, 
and  slope  of  the  applied  voltage).  We  see  that  under  these 
conditions,  the  force  is  only  slightly  dependent  on  the  dielectric 
layer,  tends  to  increase  with  the  voltage  slope  and  is  on  the  same 
order  or  even  higher  than  the  force  that  was  calculated  in  the 
case  of  a  positive  ramp  voltage  (see  figure  6). 

3.2.3.  Expansion  velocity  of  the  negative  space  charge  above 
the  surface.  The  velocity  of  the  expansion  of  the  negative 
ion  cloud  above  the  surface  can  be  defined  as  the  velocity 
of  the  relative  maximum  of  the  electric  field  of  figure  13. 
Figure  15  shows  this  velocity  as  a  function  of  the  slope  of 
the  ramp  voltage,  r]y  The  frequency  of  the  current  pulse  is 
also  plotted  as  a  function  of  r\y.  We  see  that  the  ion  cloud 
velocity  expansion  increases  linearly  with  r]V  and  is  on  the 
order  of  a  few  1 00  m  s_ 1 .  The  increase  in  the  ion  cloud  velocity 
with  Tjv  is  due  to  the  fact  that  the  charging  rate  of  the  surface 
increases  with  rjy  which  is  consistent  with  the  increase  in  the 
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Figure  15.  Frequency  of  the  negative  current  pulses  and  expansion 
velocity  of  the  ion  cloud  as  a  function  of  the  negative  voltage  slope 
for  two  values  of  the  permittivity  of  the  dielectric  layer  and  for 
simulation  domain  dimensions  L  =  8  mm,  w  =  1  mm,  h  =  3  mm, 

/  =  0.5  mm. 

current  pulse  frequency  in  figure  15.  Electrons  travel  a  larger 
distance  above  the  surface  for  successive  current  pulses  (see 
figure  11).  They  generate  negative  ions  along  this  path,  and 
this  contributes  to  the  expansion  of  the  ion  cloud  above  the 
surface.  Note  that  they  also  generate  positive  ions  along  their 
path  above  the  surface,  but  more  negative  ions  are  created  at 
the  extremity  of  the  ion  cloud  where  the  field  is  below  that  for 
which  ionization  balances  attachment. 

The  difference  between  the  negative  and  the  positive 
regime  is  that  the  expansion  of  the  negative  ion  cloud  is 
continuous  during  the  negative  regime,  whereas  the  extent  of 
the  positive  ion  cloud  is  limited  by  high  current  breakdown 
in  the  positive  regime  (i.e.  in  the  positive  regime,  the  ion 
cloud  expansion  must  restart  from  zero  after  each  current  pulse, 
which  is  not  the  case  in  the  negative  regime).  We  discussed 
in  figure  7  how  the  expansion  length  of  the  positive  ion  cloud 
depends  on  the  discharge  parameter  in  the  positive  regime. 
In  the  negative  voltage  case,  we  can  also  study  the  expansion 
length  of  the  negative  ion  cloud.  Since  this  length  is  not  limited 
by  breakdown,  we  expect  that  the  only  limitation  of  the  ion 
cloud  expansion  will  be  the  finite  duration  of  the  sinusoidal 
voltage  cycle.  In  the  conditions  of  the  discussion  above,  a  ramp 
voltage  was  assumed,  and  we  found  (see  figure  15)  that  the 
ion  velocity  Vion  increases  linearly  with  r]V  •  From  figure  15, 
we  can  write  V[on  =  a?jv,  with  a  ~  0.4mmkV-1.  We  can 
extrapolate  this  result  to  the  case  of  a  sinusoidal  voltage  by 
writing  that  the  average  slope  for  a  sinusoidal  voltage  is  given 
by  r\v  —  4FV  where  F  and  V  are  the  voltage  frequency  and 
amplitude.  Therefore  the  ion  velocity  expansion  can  be  written 
as  Vion  =  4aFV.  Assuming  that  the  negative  discharge  takes 
place  during  about  1/4  of  the  voltage  cycle  duration  (1/4 F), 
we  find  that  the  total  expansion  length  Lcioud  of  the  ion  cloud 
over  one  cycle  is  actually  independent  of  the  frequency  and 
is  given  by  Lcioud  =  aV .  With  the  numbers  corresponding  to 
figure  15  (< a  ~  0.4mmkV-1),  this  gives  Lcioud  (mm)  ~  0.4  Vkv. 
Therefore  we  find  that  the  maximum  expansion  length  of  the 
negative  ion  cloud  for  a  sinusoidal  voltage  of  amplitude  e.g. 
25  kV  is  on  the  order  of  1  cm.  These  results  are  consistent 
with  the  optical  measurements  of  the  discharge  extension  as 
a  function  of  voltage  amplitude  performed  by  Takizawa  et  al 


[13].  These  measurements  give  an  extension  length  of  about 
4  mm  at  1 0  kV  (as  in  the  model)  but  the  increase  in  the  extension 
length  with  voltage  in  the  experiments  is  faster  than  predicted 
by  the  model. 

It  is  interesting  to  note  that  in  the  model,  the  negative 
ion  cloud  expansion  length  increases  with  applied  voltage 
amplitude  in  the  negative  regime,  while  the  positive  ion  cloud 
expansion  in  the  positive  regime  decreases  with  applied  voltage 
for  a  given  frequency  and  dielectric  layer  properties  as  seen  in 
figure  7.  Note  that  in  the  optical  measurements  of  Takizawa 
et  al  [13]  the  measured  extension  length  is  the  same  in 
the  positive  and  negative  phases.  This  is  because  the  light 
detected  in  the  positive  phase  corresponds  to  the  filamentary 
discharge  and  not  to  the  position  of  the  ion  cloud.  During  the 
negative  phase,  the  measured  light  emission  corresponds  to 
electrons  spreading  along  the  surface.  As  shown  in  the  model 
calculations  the  electrons  drift  across  the  space  charge  cloud 
at  each  pulse.  Therefore,  the  measurement  of  the  extension  of 
the  light  emission  in  the  negative  phase  provides  an  estimation 
of  the  ion  cloud  length. 

4.  Conclusion 

We  have  studied  the  properties  of  a  surface  DBD  in  air  in  the 
conditions  that  are  used  in  plasma  actuators  for  flow  control. 
In  the  plasma  actuator  geometry,  one  electrode  is  below  the 
dielectric  layer,  while  the  other  electrode  is  above.  The 
discharge  takes  place  in  the  volume  above  the  surface,  between 
the  tip  of  the  top  electrode  and  the  dielectric  layer  surface. 
Linearly  increasing,  positive  and  negative  voltages  (i.e.  with 
anode  or  cathode  above  the  dielectric  layer,  respectively)  have 
been  considered  in  order  to  simplify  the  interpretation  of  the 
results.  The  complex  chemistry  of  air  has  not  been  considered 
since  the  aim  of  this  paper  was  mainly  to  get  a  better  insight 
in  the  physics  of  these  surface  discharges,  and  not  to  provide 
accurate  predictions.  The  conclusions  that  can  be  drawn  from 
this  work  are  listed  below. 

(a)  The  current  waveform  in  the  positive  discharge  (anode 
above  the  surface)  consists  of  high  current  peaks, 
associated  with  the  development  of  a  ‘streamer  like’ 
discharges  along  the  surface.  Low  current,  ‘corona 
like’  discharges  takes  place  between  the  current  pulses. 
The  current  pulse  frequency  is  on  the  order  of  a  few 
100  kHz  under  typical  conditions.  During  this  low  current 
phase,  a  positive  ion  cloud  forms  and  expands  along  the 
surface  until  breakdown  occurs.  The  expansion  length 
of  the  positive  ion  cloud  above  the  surface  is  limited  by 
breakdown  and  increases  when  the  slope  of  the  applied 
voltage  increases.  It  varies  from  1  mm  to  a  few  millimetres 
under  typical  conditions.  The  EHD  force  for  a  positive 
discharge  is  important  only  during  the  low  current  phase 
between  the  current  pulses.  The  calculated  space  and  time 
averaged  EHD  force  parallel  to  the  surface  per  unit  length 
of  the  electrode  is  on  the  order  of  O.lmNm-1  and  is 
not  very  sensitive  to  the  conditions  (for  a  long  enough 
duration  of  the  discharge).  The  results  in  air  for  the 
positive  discharge  are  similar  to  those  obtained  in  pure 
nitrogen. 
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(b)  The  negative  discharge  (cathode  above  the  surface)  in 
air  has  different  properties.  The  discharge  current  is 
composed  of  higher  frequency  (few  megahertz)  lower 
amplitude  current  pulses  (amplitude  about  10  ten  times 
smaller  than  in  the  positive  case).  A  negative  space  charge 
forms  and  expands  continuously  above  the  surface  during 
the  current  pulses  (and  its  expansion  is  not  limited  by 
breakdown).  At  each  current  pulse,  a  plasma  forms  in 
the  vicinity  of  the  electrode  tip  and  an  electron  current 
flows  above  the  surface  across  the  ion  cloud,  creating  more 
negative  ions  at  the  tip  of  the  ion  cloud,  and  charging  the 
dielectric  surface  below  and  ahead  of  the  ion  cloud.  This 
mechanism  is  responsible  for  the  negative  space  charge 
expansion  along  the  surface.  Due  to  the  presence  of  the 
negative  ion  space  charge,  an  important  EHD  force  is 
also  present  in  the  negative  discharge  above  the  surface, 
and  is  in  the  same  direction  as  the  EHD  force  generated 
by  the  positive  discharge.  The  order  of  magnitude  of 
the  EHD  force  in  the  negative  discharge  can  be  larger 
than  the  force  in  the  positive  regime,  depending  on  the 
conditions.  The  velocity  of  the  negative  cloud  expansion 
above  the  surface  increases  linearly  with  the  slope  of 
the  applied  voltage  (like  in  the  positive  case),  and  its 
value  is  on  the  order  of  a  few  100ms-1  (several  times 
smaller  than  in  the  positive  case  for  the  same  conditions). 
In  contrast  with  the  positive  case  the  expansion  of  the 
ion  cloud  is  only  limited  by  the  duration  of  the  negative 
phase.  For  a  sinusoidal  applied  voltage,  the  maximum 
expansion  length  of  the  negative  cloud  increases  linearly 
with  the  voltage  amplitude  (independently  of  frequency) 
and  is  on  the  order  of  1  cm  for  a  voltage  amplitude 
of  25  kV. 

Work  is  in  progress  in  the  following  directions: 

•  systematic  parametric  study  of  the  EHD  force  for  a 
sinusoidal  voltage  waveform, 

•  sensitivity  analysis  of  the  results  to  different  parameters 
(secondary  emission  coefficient,  ion  mobilities,  air 
chemistry), 

•  systematic  study  of  the  accuracy  of  the  numerical  scheme, 

•  coupling  of  the  plasma  model  with  a  Navier-Stokes  solver 
to  describe  the  boundary  layer  modification  and 

•  experimental  measurement  to  confirm  the  conclusions  of 
the  model. 
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We  present  a  parametric  study  of  the  electrohydrodynamic  force  generated  by  surface  dielectric 
barrier  discharge  plasma  actuators  in  air  for  sinusoidal  voltage  waveforms.  The  simulation  results 
confirm  that  momentum  is  transferred  from  the  charged  particles  to  the  neutral  species  in  the  same 
direction  during  both  positive  and  negative  parts  of  the  cycle.  The  momentum  transfer  is  due  to 
positive  ions  during  the  positive  part  of  the  cycle  (electrode  above  the  dielectric  layer  is  the  anode), 
and  to  negative  ions  during  the  negative  part  of  the  cycle.  The  relative  contribution  of  the  positive 
and  negative  parts  of  the  cycle  depends  on  the  voltage  amplitude  and  frequency.  The  model  predicts 
that  the  contribution  of  negative  ions  tends  to  be  dominant  at  low  voltage  frequencies  and  high 
voltage  amplitudes.  ©  2009  American  Institute  of  Physics.  [DOI:  10.1063/1.3183960] 


I.  INTRODUCTION 

Electric  discharges  along  the  surface  of  an  airfoil  can 
exert  significant  forces  in  the  boundary  layer  of  a  flow  (see, 
e.g.,  the  review  of  Ref.  1).  One  of  the  possible  mechanisms 
of  action  of  a  discharge  plasma  on  a  flow  is  through  the 
electrohydrodynamic  (EHD)  force  associated  with  the  mo¬ 
mentum  transfer  from  charged  particles  to  neutral  molecules, 
leading  to  the  well  known  “ion  wind.”  In  the  context  of 
boundary  layer  modification,  ion  wind  can  be  generated  in 
surface  corona  discharges  or  in  surface  dielectric  barrier  dis¬ 
charges  (surface  DBDs  or  SDBDs).  Plasma  actuators  based 
on  surface  DBDs  have  been  considered  for  aerodynamic  ap- 
plications  since  the  pioneering  work  of  Roth  et  al.  ’  Surface 
discharges  can  generate  airflow  with  velocities  less  than  10 
m/s,  while  active  airflow  control  with  SDBDs  has  been  dem¬ 
onstrated  for  low  subsonic  velocities  (up  to  30  m/s).1 

Many  papers  have  been  devoted  to  the  modeling  of  the 
EHD  force  generated  by  a  SDBD  (see,  e.g.,  Refs.  4-9).  A 
clear  understanding  of  the  conditions  optimizing  the  force 
and  of  the  maximum  EHD  force  that  can  be  obtained  with 
surface  DBDs  is,  however,  still  needed.  Although  it  has  been 
shown  by  time  resolved  measurements10-12  and  simulations 
that  both  positive  and  negative  ions  can  contribute  to  the 
total  EHD  force  in  SDBDs  in  air  and  for  an  asymmetric 
electrode  arrangement,  there  is  no  published  parametric 
study  of  the  relative  contribution  of  each  type  of  ions  to  the 
total  force,  as  a  function  of  voltage  amplitude  and  frequency. 

In  previous  papers6,7,9  we  considered  the  simplified  case 
of  linearly  increasing  voltage  waveforms  for  positive  and 
negative  discharges.  This  approach  was  useful  to  better  un¬ 
derstand  the  physics  and  to  simplify  the  interpretation  of  the 
results.  This  work  has  clearly  confirmed  that  the  EHD  force 
in  surface  DBD  conditions,  such  as  those  of  most 
experiments,  ’  is  due  to  momentum  transfer  from  ions 
to  neutral  molecules  in  unipolar  regions  of  the  discharge  (i.e., 
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regions  where  a  positive  or  negative,  non-neutral  ion  cloud 
develops  above  the  dielectric  surface).  The  momentum  trans¬ 
fer  takes  place  mainly  during  the  low  current  phases  between 
high  current  pulses  and  the  streamers  developing  along  the 
surface  do  not  contribute  significantly  to  the  total  EHD  force. 

In  this  paper  we  consider  sinusoidal  voltage  waveforms 
and  study  the  influence  of  the  voltage  amplitude  and  fre¬ 
quency  on  the  discharge  properties  and  on  the  generated 
EHD  force. 


II.  SIMULATION  OF  SDBDS  FOR  SINUSOIDAL 
WAVEFORMS 

In  this  section  we  briefly  recall  the  principles  of  the 
model  (Sec.  II  A),  then  we  show  the  evolution  of  the  current 
and  charged  particle  densities  in  a  typical  sinusoidal  SDBD 
(Sec.  II B)  of  the  EHD  force  (Sec.  II  C)  and  of  the  surface 
charge  (Sec.  II  D). 

A.  Model 

The  discharge  model  has  been  described  in  previous  pa¬ 
pers,  and  we  only  recall  here  its  main  features.  Electron  and 
ion  (positive  and  negative)  fluid  transport  equations  are 
solved  together  with  Poisson’s  equation  for  the  electric  held 
in  a  two-dimensional  (2D),  Cartesian  geometry. 

We  consider  only  one  type  of  positive  ion  and  one  type 
of  negative  ion  with  a  basic  chemistry  including  ionization, 
attachment,  and  recombination. 

The  continuity  equations  for  electrons,  positive,  and 
negative  ions  are  written,  respectively,  as 

dn  -> 

-f  +  V-  re  =  (a-  vWe\[-repnenp,  (1) 

dn  ^ 

+  V  •  =  aflrj  -  repnenp  -  rnpnnnp,  (2) 
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-^  +  V-r„=  ?7||r  e\\-rnpnnnp,  (3) 

where  ne,np,nH  are  the  electron,  positive  ion,  and  negative 
ion  densities,  Te,Tp,Tn  the  charged  particle  fluxes,  a  and  77 
the  ionization  and  attachment  coefficients  in  air,  and  rep  and 
rnp  the  electron-ion  and  ion-ion  recombination  coefficients. 

In  the  drift-diffusion  approximation,  the  charged  particle 
momentum  equations  are  equivalent  to  writing  that  the 
charged  particle  fluxes  are  the  sum  of  a  drift  term  and  a 
diffusion  term  as  in 


L 

FIG.  1.  Simulation  domain:  L= 8  mm,  w=  1  mm,  and  h= 3  mm. 


T'e  =  fiey-neE^‘^j£^ne  j,  (4) 

G  =  Vp(npE  ~  ’  (5) 

G  =  /“»(- nnE  -  ^pVn„j,  (6) 

where  (jme , jmp , jmn)  and  ( Te,Tp,Tn )  are  the  charged  particles’ 
mobility  and  temperature,  respectively,  and  kB  is  the  Boltz¬ 
mann  constant.  The  data  for  air  (ionization,  attachment,  and 
mobility  as  a  function  of  reduced  electric  field  and  recombi¬ 
nation  coefficients)  are  the  same  as  those  of  Ref.  9.  The 
charged  particle  temperatures  in  Eqs.  (4)-(6)  are  supposed  to 
be  constant  and  equal  to  1  eV  for  electrons  and  to  ambient 
temperature  for  ions. 

Equations  (l)-(6)  above  must  be  coupled  to  Poisson’s 
equation  for  the  electric  field, 

V  •(£,£)  =  —  {n  -ne-nn)  +  —  Ss,  (7) 

eo  so 

where  e0  is  the  vacuum  permittivity,  er  is  the  relative  permit¬ 
tivity  (equal  to  1  in  the  discharge  volume  above  the  dielectric 
surface  and  supposed  to  be  equal  to  5  inside  the  dielectric 
layer),  and  crSs  represents  the  contribution  of  the  charges 
deposited  by  the  discharge  on  the  dielectric  surface.  These 
charges  are  obtained  by  time  integrating  charged  particle 
fluxes  to  the  surface.  Equation  (7)  is  solved  for  the  electric 
potential  V  (E=-VV).  The  boundary  conditions  for  Poisson’s 
equation  are  zero  field  perpendicular  to  the  walls  of  the 
simulation  domain  and  potential  difference  between  the  elec¬ 
trodes  of  the  form  V=V0  cos  {cot). 

The  boundary  condition  for  the  charged  particle  densities 
on  the  walls  of  the  simulation  domain  is  zero  density  gradi¬ 
ents,  i.e.,  no  flux  (since  the  perpendicular  field  is  zero), 
which  is  equivalent  to  a  symmetry  condition. 

Finally,  the  dielectric  surface  is  supposed  to  emit  second¬ 
ary  electrons  under  ion  bombardment,  with  a  secondary 
emission  coefficient  y  equal  to  0.05, 

fe-u±  =  -  yf~p-ux,  (8) 

The  physical  model  above  is  typical  of  atmospheric  dis¬ 
charge  models  and  is  similar  to  the  model  of  Likanskii  et 
al .,4  in  the  context  of  DBD  actuators. 

All  the  results  presented  in  this  paper  have  been  obtained 
with  a  (8  X  4  mm2)  simulation  domain  shown  in  Fig.  1.  The 


grid  spacing  is  20  /zm,  i.e.,  a  uniform  mesh  of  400  X  200  has 
been  used.  In  these  conditions  the  computation  time  on  a 
personal  computer  (PC)  workstation  is  on  the  order  of  10 
CPU  hours  for  100  /zs  (this  is  only  an  order  of  magnitude, 
and  the  CPU  time  actually  depends  on  the  voltage  amplitude 
and  frequency),  i.e.,  it  takes  a  few  days  of  CPU  time  to 
simulate  several  cycles  at  a  voltage  frequency  of  5  kHz,  for 
example.  This  is  relatively  fast  due  to  the  semi-implicit  time 
integration  of  the  transport-Poisson  equations  and  to  the  fact 
that  the  integration  time  step  for  the  electron  and  ion  trans¬ 
port  equations  can  be  larger  than  the  Courant  Friedrich  Lewy 
time  (first  order  accuracy  in  time).  The  question  of  the  accu¬ 
racy  of  the  simulation  has  been  discussed  in  Ref.  6,  and 
results  from  the  model  presented  here  compare  satisfactorily 
with  results  from  a  more  accurate  numerical  model.8,19  More 
systematic  comparisons  will  be  discussed  in  details  in  a 
forthcoming  paper. 

Results  are  presented  below  for  sinusoidal  voltages 
across  the  electrodes,  with  amplitudes  and  frequencies  in  the 
ranges  of  5-30  kV  and  1-10  kHz,  respectively.  In  the  rest  of 
the  paper  we  call  positive  (or  positive  discharge)  the  part  of 
the  cycle  when  the  electrode  above  the  dielectric  surface 
plays  the  role  of  an  anode,  and  negative  (or  negative  dis¬ 
charge)  the  part  of  the  cycle  when  the  top  electrode  is  a 
cathode. 

B.  Evolution  of  a  typical  discharge 

We  present  in  this  section  the  space  and  time  evolution 
of  a  typical  SDBD  in  air  at  atmospheric  pressure.  Since  the 
time  evolution  of  a  SDBD  can  be  quite  complex,  most  of  the 
results  (time  evolution  of  the  charged  particle  densities  and 
EHD  force)  will  be  shown  as  a  function  of  time  and  position 
along  the  dielectric  surface  and  integrated  in  the  direction 
perpendicular  to  the  dielectric  surface.  This  allows  a  simpler 
visualization  of  the  time  variations  of  those  quantities  and 
provides  a  clearer  understanding  of  the  mechanisms.  This 
more  concise  representation  contains  the  important  features 
and  is  sufficient  for  our  purpose.  The  detailed  time  evolution 
of  the  2D  quantities  was  described  in  previous  papers  in  the 
case  of  linearly  increasing  voltages  (see,  e.g.,  Refs.  6-9). 

Figure  2  shows  the  time  evolution  of  the  calculated  cur¬ 
rent  for  a  voltage  waveform  of  15  kV  amplitude  and  10  kHz 
frequency.  Figures  3-5  show  the  calculated  time  variations 
of  the  electron,  positive  ion,  and  negative  ion  number  densi¬ 
ties  along  the  dielectric  surface  and  integrated  along  the  di¬ 
rection  perpendicular  to  the  surface  (therefore  given  in  units 
of  cm-2),  as  a  function  of  time,  in  the  same  conditions. 
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FIG.  2.  (Color  online)  Time  variations  of  the  calculated  current  and  applied 
voltage  for  a  15  kV  amplitude,  10  kHz  frequency  voltage  waveform.  The 
electrode  above  the  dielectric  layer  plays  the  role  of  an  anode  during  the 
“positive  phase”  (or  positive  part  of  the  cycle)  and  plays  the  role  of  a 
cathode  during  the  “  negative  phase.” 

We  see  on  the  calculated  current  in  Fig.  2  that  the  posi¬ 
tive  and  negative  parts  of  the  cycle  are  quite  different,  in 
agreement  with  the  calculations  for  positive  and  negative 
ramp  voltages  described  in  previous  papers.9  A  few  high  cur¬ 
rent  pulses  separated  by  low  current  phases  are  apparent  dur¬ 
ing  the  positive  part  of  the  cycle  (Fig.  2).  These  current 
pulses  are  separated  by  several  10  s  of  microseconds.  The 
negative  part  of  the  cycle  is  composed  of  a  much  larger 
number  of  current  pulses  (frequency  in  the  megahertz  range) 
of  smaller  amplitudes.  The  calculated  current  waveform  and 
the  asymmetry  between  positive  and  negative  parts  of  the 
cycle  predicted  by  the  model  are  in  excellent  qualitative 
agreement  with  the  experiments  (see,  e.g.,  Refs.  13-18). 

The  time  variations  of  the  charged  particle  densities 
along  the  surface  and  integrated  in  the  direction  perpendicu¬ 
lar  to  the  surface  (Fig.  3-5)  clearly  show  the  streamer  phases 
associated  with  the  high  current  pulses  in  the  positive  regime 
and  the  glow  corona  phase  during  the  negative  part  of  the 
cycle.  The  high  frequency  oscillations  in  the  negative  phase 
are  not  apparent  on  Figs.  3-5  because  the  time  resolution  of 
the  plotted  densities  is  larger  than  the  period  of  these  oscil¬ 
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FIG.  3.  (Color  online)  Time  evolution  of  the  electron  density  along  the  x 
direction  and  integrated  in  the  y  direction  (see  Fig.  1)  and  therefore  given  in 
units  of  1010  cm-2  for  the  conditions  of  Fig.  2  (10  kHz,  15  kV).  The  dashed 
line  indicates  the  time  evolution  of  the  applied  voltage  on  the  top  electrode 
(the  electrode  below  the  dielectric  layer  being  grounded). 


10‘3  10‘2  10'1  1 


FIG.  4.  (Color  online)  Time  evolution  of  the  positive  ion  density  along  the 
x  direction  and  integrated  in  the  y  direction  in  units  of  1010  cm-2  for  the 
conditions  of  Fig.  2  (10  kHz,  15  kV).  The  dashed  line  indicates  the  time 
evolution  of  the  applied  voltage  on  the  top  electrode. 

lations  (see  Ref.  9  for  more  details).  A  non-neutral,  negative 
ion  cloud  develops  along  the  surface  during  the  negative 
phase.  The  development  of  this  ion  cloud  is  consistent  with 
the  results  obtained  with  linearly  varying  negative  voltages. 
Positive  ions  are  also  formed  above  the  dielectric  surface 
during  the  negative  phase,  but  the  positive  ion  density  is 
smaller  than  the  negative  ion  density,  especially  in  the  front 
of  the  ion  cloud  (compare  Figs.  4  and  5),  and  the  space 
charge  is  negative.  The  existence  of  this  negative  space 
charge  is  possible  because  the  ion  densities  in  the  cloud  are 
small  enough  that  charge  separation  by  the  applied  field  is 
possible  in  that  region  (the  plasma  is  not  quasineutral).  The 
electric  field  in  that  region  is,  as  mentioned  in  Ref.  9,  slightly 
below  the  value  for  which  ionization  balances  attachment. 
Note  that  the  propagation  velocity  of  the  negative  ion  cloud 
along  the  surface  can  be  estimated  to  be  of  a  few  104  cm/s 
in  these  conditions. 

During  the  positive  phase  a  positive  ion  cloud  develops 
over  the  surface,  but  its  development  is  interrupted  by 
streamer  breakdown  and  restarts  from  zero  after  each 
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FIG.  5.  (Color  online)  Time  evolution  of  the  negative  ion  density  along  the 
x  direction  and  integrated  in  the  y  direction  in  units  of  1010  cm-2  for  the 
conditions  of  Fig.  2  (10  kHz,  15  kV).  The  dashed  line  indicates  the  time 
evolution  of  the  applied  voltage  on  the  top  electrode. 
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FIG.  6.  (Color  online)  Time  evolution  of  the  EHD  force  per  unit  surface  (in 
N/m2)  parallel  to  the  dielectric  surface  (force  per  unit  volume  integrated  in 
the  y  direction)  for  the  conditions  of  Fig.  2  (10  kHz,  15  kV).  The  force  is 
mainly  directed  in  the  positive  x  direction  (blue  color,  from  0  to  200  N/m2, 
contours  separated  by  20  N/m2),  i.e.,  away  from  the  top  electrode,  except 
near  the  tip  of  the  top  electrode  where  it  is  directed  toward  the  electrode  (red 
color,  maximum  of  -500  N/m2).  The  dashed  line  indicates  the  time  evolu¬ 
tion  of  the  applied  voltage  on  the  top  electrode. 

streamer  pulse.  This  is  very  different  from  the  negative  phase 
where  the  negative  ion  space  charge  is  continuous  and  is  not 
strongly  perturbed  by  the  high  frequency  current  oscillations. 

It  is  interesting  to  note  in  Figs.  3-5  that  the  length  of  the 
streamer  channel  along  the  dielectric  surface  increases  from 
the  first  to  the  third  high  current  pulses  in  the  positive  re¬ 
gime.  The  first  streamer  extends  about  3  mm  along  the  di¬ 
electric  surface,  the  second  streamer  length  is  about  6  mm, 
while  the  third  one  reaches  the  end  of  the  simulation  domain, 
i.e.,  a  length  of  7  mm  or  larger.  This  is  qualitatively  consis- 
tent  with  the  experiments  of  Allegraud  et  al.  who  measured 
the  streamer  extension  as  a  function  of  time  for  a  lower  fre¬ 
quency  (50  Hz)  SDBD  and  showed  that  the  streamer  length 
increases  with  time  during  the  voltage  rise  in  the  positive 
phase. 

C.  EHD  force 

The  instantaneous  EHD  force  per  unit  volume  is  ob¬ 
tained,  at  any  time  from  the  simulation  by5 

^ehd  =  e{ne  -np-  nn)E  -  kBTpfnp  -  kBTeVne 

-kBTn\nn.  (9) 

The  component  of  the  EHD  force  parallel  to  the  surface  and 
integrated  along  the  direction  perpendicular  to  the  surface  is 
displayed  in  Fig.  6  for  the  conditions  of  Figs.  2-5.  We  see 
that  (1)  the  EHD  force  along  the  surface  has  the  same  direc¬ 
tion  during  the  positive  and  negative  phases  of  the  sinusoidal 
cycle,  i.e.,  away  from  the  top  electrode,  (2)  there  is  a  large 
negative  force  (directed  toward  the  top  electrode)  that  is  lo¬ 
calized  in  a  very  small  region  next  to  the  top  electrode  (posi¬ 
tive  ion  sheath)  during  the  negative  part  of  the  cycle,  (3)  the 
EHD  force  during  the  negative  and  positive  phases  of  the 
cycle  have  very  different  spatial  distributions:  the  force  due 
to  negative  ions  in  the  negative  phase  is  distributed  rather 
smoothly  and  uniformly,  over  a  length  of  about  5  mm  along 
the  dielectric  surface,  while  the  force  due  to  positive  ions,  in 
the  positive  phase  of  the  cycle,  is  distributed  nonuniformly, 


negative  force 


ycuive  iuioc 

"1  T  ”1 


40 

20 


positive^  negative 
phase  \v  phase  / 


350  400  450  500 

Time  (us) 


FIG.  7.  (Color  online)  Time  integrated  total  (integrated  over  the  simulation 
domain)  force  parallel  to  the  surface,  in  the  conditions  of  Figs.  2-6  (10  kHz, 
15  kV). 

closer  to  the  top  electrode;  also  there  is  no  force  during  the 
high  current  pulses  since,  as  discussed  in  previous  papers,5,6,8 
the  EHD  force  is  quasizero  in  a  quasineutral  plasma  (the 
channel  of  the  streamers).  The  force  that  can  be  seen  during 
the  positive  phase  (seen  as  three  successive  maxima  on  Fig. 
6)  corresponds  to  the  development  of  the  positive  ion  space 
charge  before  each  streamer  pulse  (during  which  the  force 
goes  to  zero). 

The  EHD  force  of  Fig.  6  can  be  integrated  along  v  and 
up  to  time  t  and  displayed  as  a  function  of  time.  This  time 
integrated  force  is  represented  in  Fig.  7  for  the  conditions  of 
Figs.  2-6. 

We  see  in  Fig.  7  that,  in  these  conditions  (10  kHz,  15 
kV),  the  increase  in  the  EHD  force  during  the  negative  part 
of  the  cycle  is  slightly  larger  than  during  the  positive  part  of 
the  cycle.  We  find  that  the  relative  contribution  of  the  posi¬ 
tive  part  of  the  cycle  (positive  ions)  and  the  negative  part 
(negative  ions)  to  the  total  EHD  force  during  a  voltage  cycle 
strongly  depends  on  the  voltage  amplitude  and  frequency.  At 
higher  voltage  frequencies  and  lower  voltage  amplitudes,  the 
EHD  force  is  larger  during  the  positive  part  of  the  cycle 
(dominant  contribution  of  positive  ions),  while  at  lower  volt¬ 
age  frequencies  and  larger  voltage  amplitudes,  the  EHD 
force  is  larger  during  the  negative  part  of  the  cycle  (dominant 
contribution  of  negative  ions).  This  is  further  discussed  be¬ 
low. 

The  effect  of  voltage  amplitude  on  the  space  and  time 
distribution  of  the  EHD  force  for  a  given  voltage  frequency 
is  illustrated  in  Figs.  8(b)  and  9(b).  The  corresponding  volt¬ 
age  and  current  waveforms  are  shown  in  Figs.  8(a)  and  9(a). 
The  EHD  force  is  integrated  along  the  direction  perpendicu¬ 
lar  to  the  surface  and  is  displayed  as  a  function  of  position 
along  the  surface  and  time. 

Figure  8  shows  the  EHD  force  for  a  5  kHz  frequency 
and  10  kV  voltage  amplitude.  We  see  that  in  that  case,  the 
force  is  relatively  small  during  the  negative  part  of  the  cycle. 
The  spatial  extension  of  the  force  is  also  limited  to  about  3 
mm  away  from  the  top  electrode.  The  EHD  force  during  the 
positive  part  of  the  cycle  is  larger  but  its  extension  is  also 
limited  to  about  3  mm.  Note  that  under  these  conditions, 
because  of  the  low  voltage  and  relatively  low  frequency  (i.e., 
low  voltage  increase  rate),  there  is  no  streamer  formation 
during  the  positive  part  of  the  cycle,  and  the  current  during 
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FIG.  8.  (Color  online)  (a)  Calculated  current  for  a  sinusoidal  voltage  of 
amplitude  10  kY  at  a  frequency  of  5  kHz  (also  plotted),  (b)  Time  evolution 
of  the  EHD  force  per  unit  surface  (in  N/m2)  parallel  to  the  dielectric  surface 
(force  per  unit  volume  integrated  in  the  y  direction)  for  the  geometry  of  Fig. 
1  with  a  sinusoidal  voltage  of  amplitude  of  10  kV  at  a  frequency  of  5  kHz 
(a).  The  force  is  mainly  directed  in  the  positive  x  direction  (blue  color,  from 
0  to  100  N/m2,  contours  separated  by  10  N/m2),  i.e.,  away  from  the  top 
electrode,  except  near  the  tip  of  the  top  electrode  where  it  is  directed  toward 
the  electrode  (red  color,  maximum  of  -100  N/m2).  The  dashed  line  indi¬ 
cates  the  time  evolution  of  the  applied  voltage  on  the  top  electrode. 


FIG.  9.  (Color  online)  (a)  Calculated  current  for  a  sinusoidal  voltage  of 
amplitude  of  25  kV  at  a  frequency  of  5  kHz  (also  plotted),  (b)  Time  evolu¬ 
tion  of  the  EHD  force  parallel  to  the  dielectric  surface  and  integrated  in  the 
y  direction  for  the  geometry  of  Fig.  1  with  a  sinusoidal  voltage  of  amplitude 
25  kV  at  a  frequency  of  5  kHz  (a).  The  force  is  mainly  directed  in  the 
positive  x  direction  (blue  color,  from  0  to  200  N/m2,  contours  separated  by 
20  N/m2),  i.e.,  away  from  the  top  electrode,  except  near  the  tip  of  the  top 
electrode  where  it  is  directed  toward  the  electrode  (red  color,  maximum  of 
-500  N/m2).  The  dashed  line  indicates  the  time  evolution  of  the  applied 
voltage  on  the  top  electrode. 


that  phase  is  only  a  corona  current  leading  to  an  ion  cloud 
formation  above  the  surface  and  to  the  corresponding  EHD 
force. 

Figure  9  displays  the  EHD  force  for  the  same  voltage 
frequency,  5  kHz,  but  for  a  larger  voltage  amplitude,  25  kV. 
The  EHD  force  is  now  much  larger  during  the  negative  part 
of  the  cycle.  We  also  see  that  four  streamers  form  during  the 
positive  part  of  the  cycle  and  that  the  force  is  more  intense 
but  distributed  on  a  much  smaller  region  in  the  positive  part 
of  the  cycle.  The  development  of  the  EHD  force  is  clearly 
perturbed  by  the  streamers,  is  equal  to  zero  during  streamer 
formation,  and  restarts  from  zero  after  each  streamer  devel¬ 
opment. 

Note  in  Figs.  8  and  9  the  large  negative  force  (directed 
toward  the  top  electrode)  in  a  very  small  region  next  to  the 
top  electrode.  This  force  is,  as  in  Fig.  6  above,  due  to  the 
positive  ion  sheath  close  to  the  top  electrode  in  the  negative 
regime. 


D.  Surface  charge  and  potential 

The  time  variations  of  the  surface  charge  along  the  di¬ 
electric  surface  are  shown  in  Fig.  10  for  the  conditions  of 
Fig.  9  (5  kHz,  25  kV).  The  current  waveforms  are  also 
shown  on  these  figures  for  clarity. 


We  see  that  the  charge  increases  smoothly  during  the 
negative  phase.  As  we  discussed  in  previous  papers,9  the 
charge  during  the  negative  phase  is  due  mainly  to  electrons 
which  extend  further  and  further  away  from  the  top  electrode 
at  each  high  frequency  pulse.  The  high  frequency  pulses  are 
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FIG.  10.  (Color  online)  Contours  of  constant  surface  charge  density  (in  units 
of  10-7  C/m2)  on  the  dielectric  surface  as  a  function  of  time  for  the  condi¬ 
tions  of  Fig.  9  (5  kHz,  25  kV).  The  current  waveform  is  shown  below  the 
surface  charge. 
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FIG.  11.  (Color  online)  Contours  of  constant  potential  (kV)  along  the  sur¬ 
face  as  a  function  of  time  for  the  conditions  of  Fig.  9  (5  kHz,  25  kV).  The 
current  and  voltage  waveforms  are  also  shown  below  the  surface  potential. 


not  seen  in  the  figure  because  the  time  resolution  of  the  fig¬ 
ure  is  1  /x s  and  the  high  frequency  pulses  of  the  negative 
regime  are  in  the  megahertz  range.  The  charging  during  the 
positive  part  of  the  cycle  is  not  as  smooth,  and  we  see  on  Fig. 
10  that  at  each  streamer  pulse  (see  the  current  waveform 
below  the  surface  charge  plot),  positive  charges  are  abruptly 
deposited  on  the  surface  (discontinuity  in  the  slope  of  the 
surface  charge  as  a  function  of  time). 

We  also  see  on  this  figure  that  the  streamer  length  in¬ 
creases  with  time  during  a  positive  discharge. 

Figure  1 1  shows  the  time  variations  of  the  potential  dis¬ 
tribution  along  the  surface,  together  with  the  voltage  and 
current  waveforms.  The  streamers  are  also  apparent  on  the 
potential  distribution  since  the  surface  potential  below  the 
plasma  filament  generated  by  a  streamer  becomes  almost  in¬ 
stantaneously  equal  to  the  top  electrode  potential  (the  high 
current  pulses  are  very  short  on  the  time  scale  of  the  figure, 
and  the  potential  drop  along  the  filament  is  relatively  small). 


III.  PARAMETRIC  STUDY 

We  have  performed  systematic  calculations  of  the  EHD 
force  parallel  to  the  surface,  integrated  in  the  simulation  do¬ 
main,  and  averaged  in  time,  for  voltage  frequencies  between 
1  and  10  kHz,  and  voltage  amplitudes  between  4  and  30  kV. 


A.  EHD  force 

The  contributions  of  the  positive  and  negative  parts  of 
the  cycle  to  the  total  EHD  force  parallel  to  the  surface  are 
represented  in  Fig.  12  as  a  function  of  voltage  amplitude,  and 
for  different  voltage  frequencies.  The  force  is  integrated  over 
the  simulation  domain,  averaged  in  time,  and  given  per  unit 
length  of  electrode  (since  the  model  is  2D  Cartesain).  The 
total  force  parallel  to  the  surface  is  represented  in  Fig.  13  as 
a  function  of  voltage  amplitude  and  frequency. 

We  see  that,  as  mentioned  above,  the  contribution  of 
positive  ions  (positive  phase)  tends  to  be  larger  at  higher 


Voltage  Amplitude  (kV) 


FIG.  12.  (Color  online)  Calculated  EHD  force  parallel  to  the  surface,  inte¬ 
grated  over  the  simulation  domain,  and  averaged  in  time,  as  a  function  of 
voltage  amplitude,  and  for  different  voltage  frequencies.  The  contributions 
of  the  positive  (open  symbols)  and  the  negative  (full  symbols)  part  of  the 
cycle  (see  Fig.  7)  to  the  total  force  are  represented  separately.  The  EHD 
force  is  given  per  unit  length  of  the  electrode  in  the  direction  perpendicular 
to  the  simulation  domain. 


frequencies  and  lower  voltage  amplitudes  while  the  contri¬ 
bution  of  negative  ions  is  larger  at  lower  frequencies  and 
larger  voltages. 

Note,  again,  that  the  space  distributions  of  the  EHD 
forces  are  very  different  for  positive  and  negative  ions  (see 
above),  so  that  the  same  force  magnitude,  in  Fig.  12,  during 
positive  or  negative  regimes  may  not  have  the  same  aerody¬ 
namic  effect. 

The  maximum  total  force  displayed  in  Fig.  13  is  around 
100  mN/m.  Two  measurements  from  Abe  et  al .17  for  6  and 
10  kV  amplitude  voltage  sinusoidal  waveforms  at  5  kHz  are 
also  shown  for  comparisons  (taken  from  Fig.  8  of  Ref.  17, 
the  force  per  unit  length  being  deduced  by  dividing  the  mea¬ 
sured  total  force  by  the  electrode  length).  The  calculated 
force  is  about  twice  larger  than  the  measured  force.  The  cal¬ 
culated  EHD  force  is  therefore  on  the  same  order  as  the 
measured  force  and  this  can  be  considered  as  a  very  encour¬ 
aging  result  considering  the  complexity  of  the  momentum 
transfer  from  the  discharge  to  the  neutral  gas,  the  relative 
simplicity  of  the  model,  and  the  difficulty  of  the  measure- 
ments.  The  papers  by  Enloe  et  al.  and  by  Baughn  et  al. 
show  measurements  of  the  EHD  force  that  are  consistent 
with  those  of  Abe  et  al .,  i.e.,  on  the  order  of  10-20  mN/m  for 
voltage  amplitudes  around  10  kV  and  frequency  around  10 
kHz. 


Voltage  (kV) 

FIG.  13.  (Color  online)  Calculated  total  EHD  force  parallel  to  the  surface, 
integrated  over  the  simulation  domain,  and  averaged  in  time  (sum  of  the 
positive  and  negative  contributions  of  Fig.  12).  The  star  symbols  correspond 
to  the  measurements  of  Abe  et  al.  (Ref.  17)  at  5  kHz. 
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The  model  results  of  Fig.  13  also  show  that  the  force  can 
increase  significantly  when  the  applied  voltage  amplitude  is 
larger  than  10  kV.  We  shall  see  in  the  subsection  below  that 
it  is  possible  to  increase  the  voltage  to  20  or  30  kV  while 
keeping  the  dissipated  power  in  reasonable  limits,  if  one  op¬ 
erates  at  lower  voltage  frequencies. 

We  see  on  Fig.  13  a  saturation  and  a  decrease  in  the 
force  for  high  enough  applied  voltages.  This  saturation  is 
reached  at  higher  voltages  for  lower  frequencies.  The  limited 
length  of  the  simulation  domain  may  be  partly  responsible 
for  this  saturation.  Note  that  the  results  of  Figs.  12  and  13  at 
low  frequencies  (2  kHz)  and  high  voltages  should  be  re¬ 
garded  as  preliminary  since  they  were  obtained  by  simulat¬ 
ing  no  more  than  2  cycles  (about  4  days  of  CPU  time  on  a 
PC  workstation)  and  convergence  to  a  harmonic  steady  state 
can  take  several  cycles.  Nevertheless,  the  trend  predicted  by 
the  model  is  that  although  the  EHD  force  clearly  increases 
with  frequency  at  low  voltage  amplitudes,  this  is  no  longer 
the  case  at  higher  voltage  amplitudes.  This  is  to  be  compared 
with  the  experimental  results  of  Dai  and  Roth21  (see  Fig.  16c 
of  Ref.  21),  which  show  that  the  ion  wind  velocity  increases 
with  frequency  for  rms  voltages  below  about  7  kV  and  de¬ 
creases  with  frequency  above  this  voltage. 

B.  Power  dissipation 

The  averaged  power  dissipated  for  each  case  of  Figs.  12 
and  13  has  been  calculated  and  is  displayed  in  Fig.  14(a)  as 
a  function  of  voltage  amplitude  and  for  different  frequencies. 
Again,  the  results  at  high  voltage  and  low  frequency  (2  kHz) 
need  to  be  confirmed  by  calculations  over  a  larger  number  of 
cycles. 

The  EHD  force  of  Fig.  13  is  replotted  on  Fig.  14(b)  for 
comparison  with  the  dissipated  power.  Finally,  the  force  per 
unit  power  is  plotted  as  a  function  of  voltage  amplitude  and 
frequency  in  Fig.  14(c).  The  efficiency  (force  generated  for  a 
given  electric  power)  is  larger  at  lower  voltages  because 
there  are  less  high  current  pulses  (streamers)  in  these  condi¬ 
tions,  but  the  force  is  small  at  low  voltages  and  increases 
with  voltage.  It  seems  from  Fig.  14(c)  that  larger  EHD  force 
and  better  efficiency  are  obtained  at  lower  frequency  and 
high  voltage  amplitudes.  This  corresponds  to  regimes  where 
the  contribution  of  negative  ions  to  the  EHD  force  is  domi¬ 
nant.  The  best  way  to  obtain  a  large  EHD  force  with  a  rea¬ 
sonable  power  consumption  is  therefore  to  operate  at  high 
voltage  amplitude  and  low  frequency.  This  is  consistent  with 
the  experiments  of  Moreau  (see,  e.g.,  Ref.  1). 

Finally,  note  that  some  of  the  curves  presented  in  Figs. 
12-14  are  not  perfectly  smooth  (in  spite  of  the  fact  that  the 
model  is  fluid  and  there  is  no  statistical  fluctuations).  This  is 
because  the  number  of  current  pulses  in  this  2D  simulation  is 
small  (especially  in  the  positive  regime).  For  example,  when 
increasing  the  voltage  at  a  constant  frequency,  the  number  of 
streamer  pulses  during  the  positive  phase  increases  dis¬ 
cretely.  The  dissipated  power  increases  abruptly  as  a  function 
a  frequency  at  each  increase  in  the  number  of  streamers. 
However,  the  calculations  of  Fig.  14  were  not  performed 
with  small  voltage  increments  and  the  discontinuities  due  to 
the  appearance  of  new  streamers  do  not  appear  clearly.  This 


(a)  Voltage  (kV) 


FIG.  14.  (Color  online)  (a)  Power  dissipated  per  unit  length  of  the  electrode, 
(b)  total  EHD  force,  and  (c)  force  per  unit  power  as  a  function  of  voltage 
amplitudes  and  for  different  voltage  frequencies. 

explains  the  “unnatural”  change  in  slope  that  can  be  ob¬ 
served,  for  example,  on  Fig.  14(a).  This  would  not  happen  in 
a  real  experiment  because  of  the  statistical  distribution  of 
streamers  along  the  electrode. 

IV.  SUMMARY  AND  CONCLUSION 

We  have  presented  simulations  of  a  SDBD  plasma  ac¬ 
tuator  for  sinusoidal  voltages  of  frequency  and  amplitude  in 
the  1-10  kHz,  and  4-30  kV  ranges,  respectively. 

The  results  confirm  previous  calculations  performed 
with  linearly  increasing  positive  or  negative  voltages. 

The  EHD  force  is  directed  away  from  the  top  electrode 
both  in  the  positive  and  negative  parts  of  the  sinusoidal  volt¬ 
age.  The  momentum  transfer  leading  to  the  EHD  force  is  due 
to  positive  ions  during  the  positive  phase  and  to  negative 
ions  in  the  negative  phase  of  the  cycle.  The  relative  contri¬ 
butions  of  the  positive  and  negative  ions  to  the  total  EHD 
force  depend  on  the  voltage  amplitude  and  frequency.  At 
high  frequencies  and  low  voltage  amplitude,  the  contribu- 
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tions  of  positive  ions  are  dominant,  while  the  contribution  of 
negative  ions  to  the  total  EHD  force  is  dominant  at  low 
frequencies  and  high  voltage  amplitudes. 

The  space  distribution  of  the  force  is  different  for  posi¬ 
tive  and  negative  ions.  In  the  positive  regime,  the  EHD  force 
due  to  positive  ions  tends  to  be  distributed  closer  to  the  top 
electrode.  This  is  because  a  new  positive  ion  cloud  must  start 
from  the  top  electrode  after  each  streamer.  During  the  nega¬ 
tive  phase,  the  negative  ion  cloud  develops  continuously  and 
its  spatial  extension  tends  to  be  larger.  The  EHD  force  during 
the  negative  phase  is  therefore  more  spread  out  along  the 
dielectric  surface.  The  extension  of  the  force  along  the  sur¬ 
face  in  the  negative  phase  increases  when  the  frequency  de¬ 
creases  and  the  voltage  amplitude  increases.  The  differences 
in  the  space  and  time  distributions  of  the  EHD  force  during 
the  positive  and  negative  phases  may  lead  to  different  aero¬ 
dynamic  effects  even  for  the  same  total  force  (integrated  in 
space).  Finally,  a  large  negative  EHD  force  (i.e.,  directed 
toward  the  top  electrode)  distributed  in  a  small  region  at  the 
tip  of  the  top  electrode  always  exists  during  the  negative 
phase.  This  force  is  associated  with  the  positive  ion  sheath 
next  to  the  top  electrode  in  that  regime.  When  spatially  inte¬ 
grated,  this  negative  force  is  negligible  with  respect  to  the 
total  positive  force  in  the  negative  ion  cloud,  but  it  may  also 
have  some  effect  on  the  flow. 

The  force  obtained  in  the  simulation  is  on  the  order  of 
several  10  mN/m  (with  a  maximum  of  100  mN/m)  for  a 
simulation  domain  of  8  mm  length,  increases  with  voltage 
amplitude,  and  is  more  efficiently  produced  at  lower  frequen¬ 
cies. 

The  discharge  efficiency  to  produce  the  EHD  force  is 
less  than  1  mN/W.  For  a  given  frequency,  the  efficiency  in 
the  production  of  the  EHD  force  strongly  decreases  with 
voltage  amplitude  above  a  given  voltage  (because  more  en¬ 
ergy  is  spent  in  streamer  discharges  at  higher  voltage).  The 
calculated  efficiency  at  the  maximum  force  is  on  the  order  of 
0.2  mN/W,  and  this  is  consistent  with  the  measurements. 

From  the  results  presented  here  and  in  previous  papers, 
and  from  the  understanding  of  the  SDBD  provided  by  the 
model,  we  can  anticipate  what  could  be  the  effect  of  voltage 
waveforms  other  than  sinusoidal  (i.e.,  positive  or  negative 
sawtooth,13  etc.).  For  example,  we  can  expect  that  negative 
sawtooth  waveforms  at  low  frequency  and  large  voltage  am¬ 
plitudes  could  be  efficient  because  they  would  tend  to  opti¬ 
mize  the  contribution  of  negative  ions.  We  call  negative  saw¬ 
tooth  a  voltage  waveform  where  the  top  electrode  voltage 
decreases  on  a  long  time  scale  while  it  goes  back  to  zero  on 
a  short  time  scale.  However,  strong  streamers  will  of  course 
form  when  the  voltage  goes  back  to  zero  (positive  discharge) 
after  the  negative  ramp  voltage,  and  this  can  affect  the  dis¬ 
charge  efficiency  in  producing  the  EHD  force. 

The  simulation  results  presented  here  are  in  excellent 
qualitative  agreement  with  the  experiments,  i.e.,  many  of  the 
features  and  trends  observed  in  the  experiments  can  be  re¬ 
produced  by  the  model.  This  means  that  the  model  includes 
the  essential  physics,  at  least  for  the  range  of  conditions  con¬ 
sidered  in  this  paper.  One  of  the  important  results  is  the 
confirmation  that  negative  ions  are  responsible  for  the  trans¬ 
fer  of  momentum  to  the  neutral  molecules  during  the  nega¬ 


tive  part  of  the  cycle,  while  positive  ions  transfer  momentum 
during  the  positive  part.  The  asymmetry  between  the  space 
and  time  distributions  of  the  EHD  force  due  to  positive  and 
to  negative  ions,  and  the  different  relative  contributions  of 
positive  and  negative  ions  depending  on  the  voltage  ampli¬ 
tude  and  frequency  are  also  new  and  important  results  from 
the  model. 

On  a  more  quantitative  point  of  view,  the  predicted  EHD 
force  per  unit  length  and  force  per  unit  power  are  consistent 
with  experiments  although  there  is  no  published  systematic 
measurements  of  these  parameters  as  a  function  of  voltage 
amplitude  and  frequency  over  a  large  range  of  conditions 
(such  measurements  would  be  extremely  useful  for  model 
validation).  More  experimental  study  on  the  time  modulation 
of  the  EHD  force  and  ion  wind  velocity  would  also  be  ex¬ 
tremely  useful.  Published  experimental  results12,1  have 
shown  the  time  modulation  of  the  EHD  force  or  of  the  gen¬ 
erated  ion  wind,  but  no  systematic  parametric  study  with 
voltage  and  frequency  has  been  reported.  Another  point  that 
needs  clarification  is  the  apparent  saturation  of  the  EHD 
force  with  voltage  for  high  enough  voltage  amplitudes,  at  a 
given  frequency. 

The  mechanisms  leading  to  momentum  transfer  from 
charged  particles  to  neutral  species  in  surface  DBDs  are  now 
relatively  well  understood.  Further  work  needs  to  do  be  done 
to  better  understand  other  possible  effects  of  a  surface  dis¬ 
charge  on  a  flow,  i.e.,  effects  that  could  be  due  to  mecha¬ 
nisms  other  than  momentum  transfer  and  ion  wind.  Experi¬ 
ments  showing  that  filamentary  surface  discharges  generated 
by  high  voltage,  nanosecond  pulses  can  affect  a  flow  al¬ 
though  they  do  not  generate  any  ion  wind  have  been  reported 
recently.  Simulations  of  the  plasma- flow  interaction  under 
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these  conditions  will  be  presented  in  a  forthcoming  paper. 
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Abstract 

Surface  dielectric  barrier  discharges  (SDBDs)  can  modify  the  boundary  layer  of  a  flow  and  are 
studied  as  a  possible  means  to  control  the  flow  over  an  airfoil.  In  SDBDs  driven  by  sinusoidal 
voltages  in  the  l-10kHz  range,  momentum  istransferred  from  ions  to  the  neutral  gas,  as  in  a 
corona  discharge  (ion  wind),  and  the  resulting  electrohydrodynamic  force  can  generate  a  flow 
of  several  ms-1  in  the  boundary  layer  along  the  surface.  In  this  paper  we  are  interested  in  a 
different  regime  of  SDBDs  where  nanosecond  voltage  pulses  are  applied  between  the 
electrodes.  Recent  experiments  by  the  group  of  Starikovskii  have  demonstrated  that 
such  discharges  are  able  to  modify  a  flow  although  no  significant  ion  wind  can  be 
detected. 

A  two-dimensional  self-consistent  numerical  model  of  the  discharge  and  gas  dynamics  in 
conditions  similar  to  those  of  these  experiments  has  been  developed.  The  model  couples  fluid 
discharge  equations  with  compressible  Navier-Stokes  equations  including  momentum  and 
thermal  transfer  from  the  plasma  to  the  neutral  gas.  This  is  a  difficult  multi-scale  problem  and 
special  care  has  been  taken  to  accurately  solve  the  equations  over  a  large  simulation  domain 
and  at  a  relatively  low  computational  cost.  The  results  show  that  under  the  conditions  of  the 
simulated  experiments,  fast  gas  heating  takes  place  in  the  boundary  layer,  leading  to  the 
generation  of  a  'micro'  shock  wave,  in  agreement  with  the  experiments. 

(Some  figures  in  this  article  are  in  colour  only  in  the  electronic  version) 


1.  Introduction 

Significant  research  efforts  have  been  devoted  in  recent  years 
to  demonstrating  the  possibility  of  using  non-thermal  plasmas 
as  actuators  for  flow  control.  Surface  corona  discharges 
and  surface  dielectric  barrier  discharges  (SDBDs)  have  been 
extensively  studied  for  this  application.  In  the  case  of  surface 
corona  discharges  [1],  two  parallel  wires  are  placed  on  a 
dielectric  surface,  a  few  centimetres  apart.  If  a  dc  voltage 
of  several  tens  of  kilovolts  is  applied  between  the  electrodes, 
the  ion  wind  generated  by  the  corona  discharge  can  reattach 
a  flow  of  a  few  ms-1  [1],  In  a  SDBD,  one  electrode 
('top  electrode’,  or  'exposed  electrode')  is  placed  above  the 
dielectric  surface  while  the  other  electrode  is  placed  below 
the  dielectric  surface  [2-5].  In  this  asymmetric  configuration 
an  ion  wind  is  generated  from  the  tip  of  the  top  electrode, 
along  the  dielectric  surface  above  the  other  electrode.  With 


sinusoidal  voltageamplitudesfromafew  kV  to 20 or 30 kV  and 
frequencies  in  the  l-10kHz  range,  SDBDs  can  also  generate 
a  flow  of  a  few  ms-1  (up  to  10ms-1).  The  advantage  of 
SDBDs  over  surface  corona  discharges  is  that  the  dielectric 
layer  between  the  electrodes  in  the  dielectric  barrier  discharge 
configuration  limits  the  current  and  prevents  the  transition  to 
an  arc  (we  are  interested  herein  a  well-controlled,  low  power, 
non-equilibrium  discharge  regime). 

A  recent  review  of  surface  discharges  for  flow  control 
can  be  found  in  [6],  A  large  number  of  experimental 
[2-7]  and  modelling  papers  [8-15]  describing  the  plasma/flow 
interaction  in  the  sinusoidal  regime  have  been  published  in  the 
last  few  years. 

The  recent  experiments  of  Starikovskii  and  co-workers 
[16-18]  introduced  a  different  type  of  SDBD  using  a  nanosec¬ 
ond  voltage  pulse  generator  (the  voltage  pulse  can  be  several 
tens  of  kilovolts  with  rise  and  decay  time  on  the  order  of  or  less 
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than  10  ns).  U  nder  these  conditions,  the  corona  regime  that  is 
responsible  for  the  ion  wind  in  sinusoidal  regimes  (see,  e.g. 
[14, 15])  is  no  longer  present  and  the  main  discharge  regime 
is  a  streamer  regime.  This  was  confirmed  by  the  quasi-zero 
ion  wind  measured  by  Starikovskii  et  al.  However,  the  au¬ 
thors  of  [16-18]  showed  that  this  kind  of  discharge  was  able 
to  affect  the  aerodynamic  properties  of  a  flow  along  the  sur¬ 
face,  and  that  a  detached  flow  could  be  reattached  when  a 
nanosecond  voltage  pulse  was  applied  between  the  electrodes 
at  a  repetition  rate  of  a  few  kilohertz  (for  spanwise  as  well 
as  streamwise  configurations  of  the  DBD  actuators  with  re¬ 
spect  to  the  flow  direction).  One  important  conclusion  of 
the  experiments  of  Starikovskii  et  a  I  [16-18]  is  that  a  large 
part  of  the  plasma  energy  is  released  into  gas  heating  in  a 
short  time  (less  than  1  /xs),  leading  to  theformation  of  a  micro 
shock  wave. 

The  purpose  of  this  work  is  to  better  understand  and 
quantify  the  gas  dynamics  generated  by  a  nanosecond  surface 
discharge  under  conditions  similar  to  Starikovskii  et  al  using 
a  self-consistent  model  of  the  discharge  and  generated  gas 
dynamics.  A  2D  fluid  model  of  the  surface  discharge  has  been 
coupled  with  a  compressible  N  avier- Stokes  description  of  the 
gas  dynamics.  This  is  a  multi-scale  problem  in  space  and  time 
since  timescales  from  picoseconds  to  tens  of  microseconds 
and  dimensions  from  afew  micrometres  to  several  centimetres 
must  be  resolved.  A  numerical  method  specially  adapted 
to  cope  with  the  multi-scale  nature  of  this  problem  has 
been  developed.  The  method  is  based  on  asynchronous 
integration  of  the  fluid  transport  equations  [19]  (fluid  discharge 
equations  and  compressible  Navier-Stokes  equations)  using 
an  asynchronous  adaptative  mesh  refinement  (A  AM  R  method) 
technique  [20]  and  allows  to  make  a  good  compromise  between 
accuracy  and  computation  time. 

The  physical  model  is  described  in  section  2,  section  3 
presents  a  summary  of  the  numerical  approach  (which  has  been 
described  in  details  elsewhere)  and  the  simulations  results  are 
discussed  in  section  4. 

2.  Physical  model 

I  n  this  section  we  describe  the  physical  basis  of  the  discharge 
model  (section  2.1)  and  the  gas  flow  model  (section  2.2)  and 
wediscussthecoupling  between  plasma  and  flow  (section  2.3). 


2.1.  Discharge  model 

TheSDBD  model  is  based  on  a  fluid  description  of  electron  and 
ion  transport  in  air  in  the  classical  drift-diffusion  and  local  field 
approximations  that  are  generally  considered  to  be  accurate 
enough  for  discharge  models  at  atmospheric  pressure. 

Three  types  of  charged  particles  are  considered:  electrons 
(index  e  in  the  equations  below),  one  generic  type  of  positive 
ions  (index  p)  and  one  type  of  negative  ions  (index  n).  In  this 
approach  we  are  not  interested  in  the  detailed  air  chemistry 
and  taking  into  account  only  one  type  of  positive  and  negative 
ions  is  sufficient  for  our  purpose. 


The  continuity  equations  for  electrons,  positive  and 
negative  ions  are  written,  respectively,  as 


dne 

—  +  v  •  re  =  (a  -  rf) 

ot 


rep^e^p> 


(i) 


drip 


-fp 


=  a 


-  repnen p  -  r^n^n p, 


(2) 


dnn 

Ot 


Tt 


—  r  np^n^p* 


(3) 


where  ne,np,nn  are  the  charged  particle  densities,  re,  rp,  rn 
the  charged  particle  fluxes,  a,  rj,  and  rep  and  rnp  are, 
respectively,  the  ionization  coefficient,  attachment  coefficient, 
and  electron- ion  and  ion- ion  recombination  coefficients. 

I  n  the  drift-diffusion  approximation,  the  charged  particles 
momentum  equations  are  equival entto  writing  that thecharged 
particles  fluxes  are  the  sum  of  a  drift  term  and  a  diffusion  term, 
plus  a  term  due  to  a  drag  from  the  gas  flow  if  present: 


nej  +  nelt ,  (4) 

=  lip  {np^  -  +  nplt ,  (5) 

—nn~t - — +  nnTt,  (6) 

where  Oe,  mp,  fin)  and  (Te,  Tp,  Tn)  are  the  charged  particle 
mobilities  and  temperatures,  respectively  and  ks  is  the 
Boltzmann  constant;  it  is  the  velocity  of  the  gas  flow. 
The  charged  particle  mobilities  are  supposed  to  depend  on 
the  local  field  as  in  [15]  and  are  obtained  with  the  BOLSIG  + 
simulation  software  [21].  The  energy  equation  is  replaced 
by  the  local  field  approximation,  i.e.  we  assume  that  the 
energy  gained  by  the  charged  particles  from  the  electric  field 
is  locally  balanced  by  the  losses  due  to  collisions  with  neutral 
molecules.  This  rough  approximation  is  used  in  most  of  the 
atmospheric  streamer  models  [22, 23]  and  is  sufficient  for  our 
purpose.  The  ionization  and  attachment  coefficients  for  air 
in  equations  (l)-(3)  are  therefore  supposed  to  depend  on  the 
local  reduced  field  (E/n,  where«  is  the  gas  density).  Weused 
tabulated  forms  of  a  and  r]  as  a  function  of  E/n,  identical 
to  those  used  in  [15].  The  charged  particle  temperatures  in 
the  diffusion  terms  of  equations  (4)-(6)  are  taken  as  constant 
(1  eV  of  the  electrons  and  ambient  temperature  for  the  ions). 
A  ssumi  ng  an  E/n  dependence  of  the  el  ectron  temperature  may 
lead  to  un-physical  results  in  some  regions  of  the  discharge 
(e.g.  negative  glow  or  wall  sheath)  and  we  prefer  imposing 
a  constant  electron  temperature  in  equation  (4).  The  value  of 
1  eV  provides  a  reasonable  order  of  magnitude  and  we  checked 
that  the  results  presented  in  this  paper  are  not  very  sensitive  to 
this  value. 

Equations  (l)-(6)  above  must  be  coupled  to  Poisson’s 
equation  for  the  electric  field: 

•  £  if  =  e(np  —  ne  —  nn)  +  cr<5s,  (7) 


where  e  is  the  permittivity  (vacuum  permittivity,  e0,  in 
the  discharge  above  the  dielectric  surface,  or  dielectric 
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permittivity,  supposed  to  be  equal  to  5e0,  inside  the  dielectric 
layer).  <j<5s  represents thecontributionofthechargesdeposited 
by  the  discharge  on  the  dielectric  surface.  These  charges  are 
obtained  by  time  integrating  charged  particles  fluxes  to  the 
surface. 

Finally,  the  dielectric  surface  is  supposed  to  emit 
secondary  electrons  under  ion  bombardment,  with  a  secondary 
emission  coefficient  y  equal  to  0.05: 

re  •  =  -yrp  -it±, 

where  it±  is  a  unit  vector  perpendicular  to  the  dielectric 
surface.  Although  we  found  that  the  results  presented  here 
are  not  very  sensitive  to  the  value  of  the  secondary  emission 
coefficient,  a  systemati c  study  of  the  i  nfl  uence  of  the  secondary 
emission  coefficient  and  pre-ionization  density  is  beyond  the 
scope  of  this  paper  and  is  left  for  a  future  work. 


2.2.  Compressible  Navier- Stokes  equations  for  the 
neutral  gas 


The  discharge  and  gas  dynamics  equations  are  coupled 
through  different  mechanisms.  The  gas  dynamics  is  modified 
by  the  plasma  because  of  (1)  momentum  transfer  from 
charged  particles  to  neutral  species  (the  electrohydrodynamic 
(EHD)  force,  which  is  a  source  term  in  the  Navier-Stokes 
equations)  and  (2)  energy  transfer  from  charged  particles  to 
neutral  molecules  (heating  term  in  the  neutral  species  energy 
equation).  In  turn,  the  discharge  can  be  modified  by  the 
flow  because  of  (1)  the  drag  on  the  charged  particles  due  to 
the  neutral  flow  velocity  and  (2)  the  changes  in  gas  density 
possibly  associated  with  the  gas  heati  ng  and  gas  dynamics  and 
that  can  modify  the  charged  particle  transport  coefficients  and 
ionization  rate. 

These  effects  are  all  included  in  the  developed  model. 
In  the  conditions  of  the  simulations  presented  in  this  paper, 
the  plasma  generates  some  gas  dynamics  effects,  but  the  gas 
dynamics  do  not  affect  the  discharge  because  a  very  short 
discharge  pulse  is  considered  (i.e.  the  gas  density  does  not 
change  significantly  during  the  discharge  even  though  gas 
heating  can  be  large).  On  the  other  hand,  no  external  flow  is 
taken  into  account  in  the  results  presented  here.  The  coupling 
between  the  discharge  model  and  the  gas  dynamics  model 
is  therefore  mainly  due  to  the  EHD  force  per  unit  volume 
Fehd;  which  can  impart  momentum  to  the  neutral  species  and 
generate  a  gas  flow,  and  to  the  power  density  Pt h  associated 
with  gas  heating  by  the  charged  particles,  which  can  generate 
pressure  gradients  and  thus  gas  flow.  The  forms  of  FEHd  and 
Pth  arediscussed  in  section  2.3.  Wefirstwritethecompressible 
N  avi  er-  Stokes  equati  ons  for  the  neutral  gas  i  ncl  udi  ng  the  E  H  D 
term  and  the  heating  term. 

The  continuity  equation  for  the  neutral  species  can  be 
written  as 


dp 

dt 


+  lt-(plt)  =  S  =  M 


+repMe»p+Pip«nrapJ, 

(8) 


where  p  is  the  mass  per  unit  volume  of  the  neutral  gas,  M  the 
mass  of  an  air  molecule  [p  =  Mn,  where  n  is  the  gas  density) 
and  it  is  the  average  velocity  of  the  flow. 


The  momentum  equation  can  be  written  as 

d-^^~ .(pit  ®lt +pl -f)  =  Spit +Fhw,  (9) 

where  PEhd  is  the  EHD  force  due  to  the  discharge  (see 
section  2.3),  is  the  scalar  gas  pressure,  I  the  unity  tensor  and 
f  the  shear  stress  tensor,  whose  components  have  the  classical 
form 


fl 


2 

( djUj  +  dtUj )  -  -  -1 f )  ■ 


(10) 


where  p  is  the  viscosity  and  <5,7  the  K  ronecker  symbol. 

Finally,  the  energy  equation  for  the  neutral  species  is 
written  as 

^  •  \(pE  +  p)  it  -  k -  f  •  it] 
dt  L  J 

=  SpE  +  PEhd  •  it  +  Pth,  (11) 

where  E  is  the  total  specific  energy  (E  =  (1  /(y  -  1  ))kBT  + 
\u2),  T  the  gas  temperature,  k  the  thermal  conductivity  of  air, 
y  the  heat  capacity  ratio  and  P*  the  power  density  associated 
with  gas  heating  from  the  plasma  (see  section  2.3). 

Finally,  the  system  of  equations  above  is  closed  with  the 
equation  of  state  of  perfect  gases: 

P  =  "T7^b  T .  (12) 

M 

Thediff erent  constant  parameters  (thetemperaturedependence 
of  viscosity  and  thermal  conductivity  is  not  considered  here) 
in  the  above  equations  are  taken  as  p  =  1.79  x  10-5  Pas, 
k  =  0.026  W  m-1  K-1,  y  =  1.4. 


2.3.  Plasma/flow  coupling  terms 

As  mentioned  above,  the  fluid  flow  dynamics  is  perturbed  by 
the  plasma  through  the  EHD  force  term  PEHd  and  the  heating 
power  density  Pt h. 

The  EHD  force  term  has  been  discussed  in  a  previous 
paper  [11]  and  is  written  as 

Pehd  =e(ne—np  —  nn )  -  ke  Tp~^  np  -  ks  Te~^  ne 

-keTn^nn.  (13) 

The  energy  transfer  from  the  plasma  to  the  neutral  gas  is  due 
to  collisions  between  charged  particles  and  neutral  species. 
We  assume  that  the  energy  gain  by  ions  from  the  electric 
field  is  totally  and  instantaneously  transferred  into  gas  heating. 
Therefore  the  contribution  of  ion- neutral  molecule  collisions 
to  the  gas  heating  term  Pt h  is  simply  equal  to  (7p  +7n)  •  ■ 

The  thermal  energy  transfer  from  electrons  to  neutral 
species  is  more  complicated  and  this  question  has  been 
discussed  in  recent  papers  [24-26].  The  energy  transferred 
by  electrons  per  unit  volume  and  per  unit  time  to  the  gas 
molecules  through  elastic  collisions  and  rotational  excitation 
is  instantaneously  converted  into  gas  heating  and  represents 
a  small  part  (a  few  per  cent  at  low  electric  fields  [28])  of  the 
total  power  density  absorbed  by  the  electrons  Peiec  =  ere E 2- 
where  ere  is  the  electron  conductivity.  The  energy  transferred 
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by  electrons  into  other  degrees  of  freedom  of  the  molecules 
(vibrational  excitation  or  electronic  excitation)  is  released  into 
gas  heating  on  different  timescales  [28].  In  this  paper  we 
are  interested  in  pressure  waves  that  can  be  generated  by  fast 
energy  deposition  in  the  gas.  Since  a  detailed  model  of  air 
chemistry  is  outside  the  scope  of  this  paper  we  used  a  simple, 
phenomenological  model  of  energy  transfer.  It  seems  that 
there  is  an  agreement  between  different  authors  [24-26]  that  a 
relatively  large  percentage  of  the  energy  spent  by  the  electrons 
in  electronic  excitation  of  N2  is  transferred  very  quickly  into 
gas  heating  due  to  the  quenching  of  nitrogen  excited  states  by 
oxygen  molecules.  We  assume  here  that  30%  of  the  electron 
energy  spent  in  electronic  excitation  of  air  is  instantaneously 
released  into  gas  heating,  which  is  consistent  with  [28].  The 
energy  stored  in  vibrational  excitation  is  released  into  gas 
heating  through  vibration-translation  (V-T)  collisions  on  a 
longer  timescale.  Since  our  model  does  not  allow  a  detailed 
description  of  this  process,  we  use  the  energy  relaxation  time 
for  V-T  collisions,  tVt,  as  a  parameter. 

The  power  density  spent  in  gas  heating  is  therefore 
given  by 

Ah  =  (  Ai-r  +  A  +  At)  +  P\  ons-  (14) 

T  he  terms  in  parentheses  on  the  rhs  of  equation  (14)  correspond 
to  the  contribution  of  electrons  to  gas  heating  while  Ans  is 
the  contribution  of  the  ions.  The  electron  contribution  to  gas 
heating  is  divided  into  Ai-r  (elastic  and  rotational  excitation 
collisions,  supposed  to  be  instantaneously  released  into  gas 
heating),  A  (electronic  excitation,  i.e.  30%  of  the  energy  going 
in  electronic  excitation  of  air,  supposed  to  be  instantaneously 
released  into  gas  heating),  and  At  (vibrational  excitation, 
released  into  gas  heating  with  a  time  constant  tvt). 

Ai-r  and  A  are  written  as 

Al-R  =  ^el-RAlec.  Pe  =  £)?EAlec>  with  Peiec  =  CeP2, 

(15) 

where  ^i-r  is  the  fractional  power  deposited  in  elastic 
collision  and  rotational  excitation  of  air  molecules,  rj E  the 
fractional  power  deposited  in  electronic  excitation  of  the 
nitrogen  molecules  and  §  =  30%. 

At  in  equation  (14)  represents  the  power  density 
corresponding  to  the  conversion  of  vibrational  energy  into 
translational  energy  (V-T  energy  transfer)  and  is  obtained  as  a 
solution  of  the  phenomenological  equation 

3  At  1  1 

— —  + — At  =  — w  Peiec.  (16) 

ot  tvt  tvt 

where  ?yv  is  the  fractional  electron  power  deposited  in 
vibrational  excitation  of  the  air  molecules. 

Equation  (16)  above  simply  means  that  the  energy  stored 
in  vibrational  excitation  is  released  into  heating  with  a  time 
delay  equal  to  tVt- 

?7ei — r  ,  ue  and  r)M  have  been  calculated  as  a  function  of  E/n 
using  the  BOLSIG  +  software  [21]  and  are  plotted  in  figure  1. 

T  he  gas  heati  ng  model  descri  bed  above  i  s  very  si  mpl  e  and 
phenomenological  but  is  sufficient  for  our  purpose.  A  similar 
descri  ption  of  the  release  of  the  energy  stored  i  n  the  vi  brational 
manifold  into  gas  heating  has  been  used  by  Naidis  [27] 


Figure  1.  Fractional  power  dissipated  by  electrons  in  air  as  a 
function  of  E/n,  into  elastic  collisions  and  rotational  excitation, 
?7ei-R ,  electronic  excitation,  >?E,  vibrational  excitation,  ^v, 
ionization,  r?|. 

and  Benilov  and  Naidis  [28]  (except  that  we  did  not  try 
to  use  a  self-consistent,  temperature  dependent  value  of  the 
vibrational  relaxation  time  tVt,  but  we  rather  considered  tVt 
as  a  parameter). 

3.  Numerical  approach 

We  describe  here  only  the  main  features  of  our  numerical 
approach.  The  details  concerning  the  numerical  techniques 
and  discretization  schemes  that  we  used  can  be  found  in  the 
references. 

3.1.  A  multi-scale  problem 

Solving  the  equations  described  in  the  previous  section  in  real 
conditions  is  a  challenge  because  the  physical  problem  is  multi¬ 
scale  in  nature,  both  in  time  and  space.  TheSDBD  leads  to  the 
formation  of  surface  streamers  [11-15]  with  plasma  densities 
up  to  1020  irr3  or  larger.  The  M  axwell  relaxation  time,  rM  = 
(eo/cre),  is  smaller  than  10— 12  s  for  these  plasma  densities, 
which  means  that  for  explicit  time  integration  of  the  Poisson- 
transport  equations,  the  integration  time  step  dr  must  be  less 
than  tm  (explicit  time  integration  of  the  Poisson-transport 
system  of  equations  means  that  Poisson's  equation  and  charged 
particle  transport  equations  are  solved  successively  and  that  the 
electric  field  calculated  at  time  tk  is  supposed  to  be  constant 
during  the  integration  of  thecharged  particle  transport  equation 
between  times  t k  and  tk+1  =  tk+dt).  Another  constraint  if  the 
continuity  equation  issolved  explicitly  (i.e.  density  at  time  A+1 
is  expressed  explicitly  as  a  function  of  density  at  time  tk)  is 
that  the  integration  time  step  of  the  electron  transport  equations 
must  be  less  than  the  Courant- Friedrich- Levy  (CFL)  time, 
rCFL  =  (Ax/i>e)min  (Ax  is  the  grid  size  and  ve  the  electron 
mean  velocity),  which  is  also  smaller  than  10_12s  in  our 
conditions.  The  gas  motion  involves  much  longer  timescales. 

We  want  to  study  the  phenomena  over  tens  or  hundreds  of 
microseconds,  i.e.  more  than  107  times  the  smaller  time  step. 
The  grid  spacing  must  be  small  enough  to  resolve  the  space 
charge  sheath  at  the  head  of  the  streamer  or  next  to  the  exposed 
electrode  tip  when  this  electrode  plays  the  role  of  a  cathode. 
The  sheath  can  be  smaller  than  10 ^m  under  atmospheric 
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pressure  conditions.  On  the  other  hand,  the  streamer  channel 
length  can  reach  1  cm  or  more,  and  we  want  to  follow  the  gas 
perturbation  induced  by  the  discharge  over  lengths  of  several 
centimetres,  i.e.  more  than  1000  times  the  sheath  length. 

The  multi  -seal  e  nature  of  thi  s  probl  em  i  n  space  and  thefact 
that  only  small  regions  (e.g.  streamer  head  or  cathode  sheath) 
need  a  very  fine  grid  impose  the  use  of  an  adaptative  mesh 
refinement  technique. 

TheCFL  and  Maxwell  relaxation  time  constraints  could 
be  overcome  by  solving  implicitly  or  semi-implicitly  charged 
particle  transport  equations  and  the  plasma-field  coupling 
(Poisson's  equation). 

In  previous  papers  [11-15]  we  performed  implicit  time 
integration  of  the  charged  particle  transport  equations  using  the 
Scharfetter-Gummel  discretization  scheme,  associated  with 
a  semi-implicit  integration  of  Poisson's  equation  (i.e.  when 
solving  Poisson's  equation  at  time  tk,  the  charged  particle 
densities  on  the  right  hand  side  of  Poisson's  equation  are 
estimated  attirne^1).  This  allows  integration  timesteps  much 
larger  than  theCFL  and  M  axwell  relaxation  times.  However, 
the  accuracy  of  this  approach  is  sufficientonly  for  small  enough 
grid  spacing.  For  large  grid  spacing  errors  due  to  numerical 
diffusion  become  important. 

In  order  to  get  a  good  compromise  between  accuracy 
and  computation  time  on  large  simulation  domains  we  have 
devel oped  an  asynchronous  ti  me  i  ntegrati on  method  associ ated 
with  an  adaptative  mesh  refinement  technique,  the  AAMR 
method  described  below. 

3.2.  Strong  plasma/aerodynamic  coupling  using  AAMR 

As  said  above  and  because  of  the  multi-scale  nature  of  the 
problem,  it  is  necessary  to  use  an  adaptative  mesh  refinement 
technique  to  allow  large  grid  spacing  in  regions  of  the 
simulation  domain  that  are  less  active.  To  prevent  excessive 
numerical  diffusion  due  to  large  grid  spacing  the  transport 
equations  must  be  solved  explicitly  and  with  a  second  order 
scheme.  O  n  the  other  hand,  the  C  F  L  ti  me  constrai  nt  i  mposes  a 
time  integration  step  At  smaller  than  the  minimum  CFL  time 
rCFL  =  (Ax/ue)min.  which  can  be  extremely  small. 

The  local  CFL  time,  tcfl, local  =  (Ax/ve),  can  be  very 
different  in  different  regions  of  the  simulation  domain.  For 
example,  the  local  CFL  time  is  very  small  in  the  streamer 
head  where  the  grid  spacing  must  be  on  the  order  of  a  few 
micrometres  and  the  electron  velocity  can  be  larger  than 
106ms_1,  but  the  local  CFL  time  is  much  larger  in  the  less 
active  region  of  the  discharge.  In  order  to  save  computation 
time,  it  is  therefore  tempting  to  perform  time  integration  of 
the  transport  equations  with  different  time  steps  at  different 
locations  of  the  grid.  This  is  the  principle  of  the  asynchronous 
method  described  by  Unfer  etal  [19]. 

The  numerical  scheme  used  for  the  integration  of  the 
charged  particle  transport  equations  is  the  classical  M  USCL 
scheme  [29].  The  MU5CL  method  is  based  on  a  linear 
piecewise  reconstruction  of  the  interface  states  based  on  the 
cell-centred  values  giving  second  order  space  accuracy.  To 
avoid  oscillations  the  reconstructed  slopes  arelimited  by  using 
a  minmod  limiter.  The  constraint  on  the  time  step  for  this 
numerical  scheme  is  actually  half  the  CFL  value  above. 


As  shown  in  [19],  a  factor  of  10  (or  more,  depending  on 
the  geometry)  in  computation  time  can  be  gained  when  using 
asynchronous  time  integration  of  the  discharge  equations  with 
a  M  USCL  scheme,  compared  with  a  standard  M  USCL  scheme 
with  synchronous  time  integration. 

The  method  has  also  been  extended  to  solve  the  2D 
laminar  compressible  unsteady  Navier- Stokes  equations.  The 
asynchronous  integration  formalism  that  we  have  developed 
allows  to  use  directly  the  HLLE  Riemann  solver  [30]  with  a 
M  USCL  scheme  for  the  reconstruction  of  the  interface  states 
and  a  classical  viscous  flux  discretization. 

The  asynchronous  time  integration  method  has  been 
recently  extended  to  cope  with  an  adaptative  mesh  refinement 
technique,  and  we  call  this  method  AAMR  [20],  Thanks  to 
the  AAMR  technique,  strongly  coupled  calculations  can  be 
performed  using  the  same  mesh  for  plasma  and  aerodynamics. 
This  enables  to  control  without  any  interpolation  the 
momentum  and  energy  transfer  from  the  plasma  to  the  neutral 
gas.  Different  classes  of  time  steps  are  used,  for  electrons, 
ions  and  neutral  species.  The  electron  time  steps  are  also 
used  to  compute  the  source  terms  of  the  seven  conservation 
equations.  From  a  Navier- Stokes  point  of  view  this  means 
that  momentum  and  energy  can  be  accumulated  into  a  cell 
without  being  transported  immediately. 

The  adaptive  meshing  relies  on  two  criteria.  The  first 
criterion  is  based  on  a  pseudo-Debye  length  for  the  plasma: 
<5  =  /  ~  (see  [20]).  The  second  criterion  relies  on 

y  e  L nernp-t-nn) 

velocity  jumps  for  the  Navier-Stokes  equations.  In  practice 
if  any  component  of  the  neutral  gas  velocity  differs  by  more 
than  lmms-1,  the  cell  is  refined.  A  cell  is  merged  if  the 
velocity  components  differ  less  than  0.1  mm  s_1.  If  any  of  the 
refinement  criteria  are  met,  the  cell  is  refined  and  the  cells  are 
merged  only  when  both  fusion  criteria  are  met.  I  n  order  to  keep 
the  total  cell  number  under  control,  different  minimal  cell  size 
are  usedfor  each  criterion.  Theminimal  cell  sizefortheplasma 
is  smaller  than  for  the  Navier-Stokes  equations.  However,  the 
Navier-Stokes  fusion  criterion  has  to  be  met  even  for  the  cells 
which  are  smaller  than  the  minimal  Navier-Stokes  cell  size. 
This  is  mandatory  to  prevent  that  the  gradients  created  by  the 
plasma  in  the  Navier-Stokes  equations  become  averaged  out 
when  the  plasma  decays. 

The  combination  of  the  asynchronous  time  integration 
and  adaptative  mesh  refinement  (AAMR)  leads  to  much  faster 
calculations  (a  factor  of  100  or  more).  The  algorithms  and  the 
coding  must,  however,  be  carefully  optimized  to  obtain  full 
benefit  from  the  AAM  R  method. 

4.  Simulation  results 

4.1.  Computational  domain  and  initial  mesh 

The  dimensions  of  the  simulation  domain  used  in  the  rest  of 
this  paper  are  0.0768  x  0.0768  m  (see figure  2).  The  dielectric 
layer  is  300 /nm  thick.  The  upper  electrode  starts  at  0.03  m 
and  ends  at  0.0351m  (5.1mm  long).  The  lower  electrode 
starts  at  0.0351  m  and  ends  at  0.0399  m  (4.8  mm  long).  The 
electrode  thickness  is  37.5  p,m.  U  nless  otherwise  mentioned, 
theminimal  cell  size  used  is  18.75  ptm  forthe  plasma  criterion 
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(a)  x  (m)  (b)  x  (m) 

Figure  2.  Simulation  domain  and  initial  mesh:  (a)  whole  domain  (the  fine  mesh  around  the  tip  of  the  top  electrode  is  not  shown),  (b)  zoom 
of  the  region  close  to  the  tip  of  the  top  electrode,  where  the  discharge  is  initiated. 


Figure  3.  Applied  voltage  pulse  and  calculated  current  for  three 
different  values  of  the  minimum  cell  size. 

and0.15  mmfortheNavier-Stokescriterion.  Smaller  minimal 
cell  size  for  the  plasma  has  also  been  used  for  comparison. 
Figure  2  shows  the  initial  mesh,  with  small  cell  sizes  in  the 
large  electric  field  region  next  to  the  tip  of  the  electrode  above 
the  dielectric  layer. 

The  initial  charged  particledensities  (electron  and  positive 
ions)  are  supposed  to  be  uniform  and  equal  to  1010  m-3  in  the 
whole  simulation  domain. 

The  results  presented  below  have  been  obtained  for  an 
applied  voltage  pulse  of  total  duration  35  ns,  with  linear  rise 
and  decay  times  (7  ns  and  15  ns,  respectively)  and  a  14 kV 
plateau  of  13  ns.  The  applied  voltage  pulse  is  shown  in  figure3, 
together  with  the  calculated  current. 

4.2.  Plasma  dynamics  for  a  35  ns  voltage  pulse  at!4kV 

The  calculated  current  of  figure  3  is  composed  of  two 
pulses.  The  first  current  pulse  starts  during  the  voltage  rise, 
passes  through  a  first  maximum  at  about  5  ns,  slightly  decays 
(due  to  the  charging  of  the  surface),  and  increases  again 
(because  of  the  continuous  increase  in  the  applied  voltage) 
and  reaches  a  second  maximum,  on  the  order  of  103A  m_1, 
at  the  end  of  the  voltage  rise.  During  the  voltage  plateau, 
the  current  decays  because  of  the  charging  of  the  dielectric 


surface.  During  the  first  current  pulse  the  top  electrode  plays 
the  role  of  an  anode  and  a  positive  ion  sheath  (streamer  head) 
develop  along  the  surface,  charging  positively  the  dielectric 
layer  surface.  When  the  applied  voltage  starts  to  decay,  the 
current  decreases  faster,  changes  sign  and  a  second  (negative) 
current  pulse  occurs.  The  second  current  pulse  is  due  to  a 
breakdown  between  the  top  electrode,  whose  voltage  rapidly 
decreases  to  zero  (in  15  ns),  and  the  dielectric  surface,  which 
has  been  positively  charged  by  thefirst  pulse.  Even  though  the 
top  electrode  voltage  does  not  become  negative,  its  potential 
quickly  becomes  I  ess  than  the  potenti  al  of  the  di  el  ectri  c  surface 
and  a  negative  current  is  collected.  We  will  see  below  that 
the  plasma  due  to  the  first  pulse  is  still  present  when  the 
second  current  pulse  takes  place.  The  second  current  pulse 
is  therefore  due  to  collection  of  charges  generated  by  the 
previous  pulse  pi  us  charges  created  by  thestrong  multiplication 
associated  with  the  increasing  potential  drop  between  the 
dielectric  surface  and  the  top  electrode  during  the  decay  of  the 
applied  voltage.  The  top  electrode  plays  the  role  of  a  cathode 
during  this  phase.  Notethatthechargecollected  during  thefirst 
current  pulse  (integral  of  the  current  between  0  and  ~20ns) 
is  larger  than  the  charge  collected  during  the  second  pulse, 
i.e.  the  total  charge  on  the  dielectric  surface  at  the  end  of  the 
discharge  is  still  positive.  This  means  that  if  a  sequence  of 
similar  voltage  pulses  is  applied  between  the  electrodes,  the 
discharge  pulses  will  evolve  so  that  the  charges  transferred 
by  each  successive  positive  and  negative  current  pulses  are 
identical  and  the  average  charge  on  the  dielectric  surface  is 
positive. 

The  current  displayed  in  figure  3  is  plotted  for  three 
different  values  of  the  minimum  cell  size  (19,  9  and  5 /am). 
A  Ithough  slight  differences  can  be  observed,  the  agreement  is 
satisfactory. 

Figures  4  and  5  show  the  time  evolution  of  the  electron 
density  and  potential  at  three  times  of  the  first  and  second 
current  pulses,  respectively. 

During  the  first  pulse  (figure  4 (a)),  the  plasma  filament 
(streamer)  propagates  along  the  dielectric  surface  with  a 
velocity  on  the  order  of  5xl07cms_1  (3mm  in  6ns).  This 
can  be  compared  with  the  experimental  results  of  Stari kovski i 
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Figured  (a)  Electron  density  (grey  scale)  and  potential  distribution  (contours)  at  three  different  times  of  the  first  current  pulse.  The  plotted 
contour  levels  are,  in  units  of  1021  rrr3,  0.1,  0.5, 1  for  t  =  4  ns  (maximum  2  x  1021  rrr3),  0.01,  0.02,  0.05,  0.1,  0.2, ...  for  t  =  7  ns 
(maximum  1.1  x  1021  nr3)  and  for  t  =  10  ns  (maximum  0.7  x  1021  nr3).  The  electric  potential  contour  levels  are  indicated  in  kV  in  the 
figure,  (b)  Space  charge  density  distribution  (n\  -  «e  -«n)  and  potential  contours  in  the  region  of  the  streamer  sheath  at  timer  =  7  ns  under 
the  conditions  of  (a).  The  represented  contours  of  the  space  charge  density  are  2, 4, 6,  8, 10, 12  x  1019  nr3,  (c)  Computational  mesh  at 
r  =  7  ns. 


etal  [18]  who  measured  a  velocity  of  108  cm  s_1  under  similar 
conditions.  This  shows  that  realistic  simulation  results  can 
be  obtained  even  though  photoionization  or  photoemission  are 
not  included  in  the  model.  The  fast  propagation  of  the  surface 
streamer  in  the  simulation  is  due  to  the  presence  of  initial 
charges  in  the  discharge  volume  and  to  secondary  electron 
emission  by  ion  impact  (ions  generated  in  the  streamer  head 
close  to  the  surface  can  cross  a  10  /cm  sheath  on  a  nanosecond 
timescale).  The  thickness  of  the  streamer  channel  is  on  the 
order  of  100  ^m. 

The  potential  along  the  dielectric  surface  increases  during 
the  first  pulse  because  of  the  charging  of  the  dielectric  by 
positive  ions.  At  t  =  10ns  (figure  4(a)),  the  average  field 
E  along  the  surface  is  on  the  order  of  28  kV  cm-1  (14  kV 
over  5  mm),  and  E/p  is  about  37  V  cm_1Torr_1;  i.e.  slightly 
over  the  critical  field  for  which  ionization  and  attachment 
frequencies  are  equal,  in  air.  As  soon  as  the  applied  voltage 
stops  increasing,  the  field  over  the  dielectric  surface  drops 


below  the  critical  field  because  of  the  charging  of  the  surface 
by  positive  ions  and  the  plasmas  decays  (current  decay  after 
t  =  10  ns  in  figure  3). 

The  streamer  head  cannot  be  clearly  seen  in  figure  4(a) 
because  of  its  small  size.  Figure  4 (b)  shows  a  close  view  of 
the  streamer  head  region  and  displays  the  net  charge  density 
distribution  [n\  -  ne  -  nn)  and  potential  contours  at  time 
t  =  7  ns.  We  can  see  the  large  potential  drop  over  a  few  tens  of 
micrometres  at  the  tip  of  the  filament.  This  region  corresponds 
to  the  streamer  head  and  the  net  charge  density  (positive  ion 
sheath)  in  that  region  is  on  the  order  of  102°  rrr3. 

Figure  4(c)  illustrates  the  adaptive  mesh  refinement 
method  and  displays  the  mesh  used  in  the  simulation  at  time 
t  =  1  us. 

During  the  decrease  in  the  applied  voltage  (i.e.  after 
t  =  20  ns),  the  potential  on  part  of  the  dielectric  surface 
becomes  higher  than  the  top  electrode  voltage.  This  can  be 
seen  in  figure  5  where  a  maximum  of  the  dielectric  surface 
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Figure 5.  Electron  density  (grey  scale)  and  potential  distribution 
(contours)  at  three  different  times  of  the  first  current  pulse.  The 
plotted  contour  levels  are  0.01,  0.02,  0.05,  0.1,  0.2, . . . ,  xlO20  m-3, 
with  a  maximum  value  of  1.2  x  1020  nr 3  for  t  =  30  ns, 

2.4  x  1020  m-3  for  t  =  33  ns  and  2.2  x  1020  m-3  for  t  =  36  ns.  The 
electric  potential  contour  levels  are  indicated  in  kV  in  the  figure. 

potential  appears  at  a  position  x  ~  3.7cm.  Because  of  this 
potential  distribution  with  a  potential  maximum  at*  ~  3.7  cm, 
the  electron  current  during  the  second  pulse  flows  to  this  point 
on  the  dielectric  surface  from  both  sides,  i.e.  from  the  left  side 
(current  flowing  from  the  top  electrode  which  plays  the  role  of 
a  cathode)  and  from  the  right  side  (where  a  smaller,  residual 
plasma  density  from  the  first  pulse  is  still  present). 

We  see  in  figure  5  the  propagation  of  an  intense  filament 
along  the  dielectric  surface  from  the  top  electrode  (cathode) 
to  the  point  of  maximum  electric  potential  at  x  ~  3.7  mm. 
The  positive  ion  density  (not  shown  here)  shows  the  presence 
of  an  ion  sheath  on  the  tip  of  the  top  electrode.  As  said 
above,  electrons  also  flow  to  the  point  of  maximum  potential 
atx  ~  3.7mm,fromtherightsideofthispoint.  Thedielectric 
surface  in  this  region  also  plays  the  role  of  a  cathode  and  we 
can  see  a  depletion  of  the  electron  density  close  to  the  surface, 
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Figure  6.  Spatial  distribution  of  the  gas  pressure  perturbation  A p 
(with  respect  to  atmospheric  pressure)  in  the  micro  shock  wave 
generated  by  the  discharge  of  figures  3-5  at  four  different  times.  The 
pressure  range  of  the  plot  is  indicated  in  each  figure.  The  contour 
line  (AP  =  0)  separates  the  region  of  positive  and  negative  A p. 

Gas  heating  from  elastic  collisions,  rotational  to  translational  energy 
transfer  and  electronic  excitation  (30%  of  electronic  excitation 
instantaneously  transferred  into  translational  energy)  is  included. 


due  to  the  presence  of  an  ion  sheath.  The  filament  on  the  right 
side  of  the  potential  maximum  is  much  less  intense  than  the 
one  on  the  top  electrode  side. 

4.3.  G  as  dynamics  and  shock  wave  initiated  by  the  discharge 

The  Navier-Stokes  equations  are  solved  together  with  the 
plasma  equations,  during  and  after  the  discharge  pulses.  The 
effect  of  the  discharge  on  the  gas  dynamics  can  be  seen  on 
the  plots  of  the  pressure  variations  induced  by  the  discharge 
at  different  times.  Figure  6  shows  the  space  distribution  of 
the  pressure  perturbation  A P  induced  by  the  discharge  at 
four  different  times  (4,  8,  16  and  25yts).  The  perturbation 
propagates  like  a  micro  shock  wave  generated  in  a  region 
around  the  tip  of  the  top  electrode.  We  see  that  A P  is  about 
1000  Pa  when  the  perturbation  front  is  5  mm  away  from  the 
discharge  region  where  it  has  been  generated.  The  results  of 
figure  6  have  been  obtained  assuming  that  the  contributions 
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Figure 7.  Spatial  distribution  of  the  gas  pressure  perturbation  at 
timer  =  8/xs,  in  the  micro  shock  wave  generated  by  the  discharge 
under  the  conditions  of  figures  3-5,  and  with  different  assumptions 
on  gas  heating:  [a)  no  gas  heating,  (b)  only  heating  by  ions  is 
included,  (c)  ions,  electron  elastic  and  rotational  excitation  and 
electronic  excitation  (30%  instantaneously  converted  into  gas 
heating)  and  (d)  same  as  (c)  with  vibrational  excitation  included, 
with  a  relaxation  time  for  V-T  collisions  tVt  =  1/xs.  The  pressure 
range  of  the  plot  is  indicated  in  each  figure. 


±10  Pa  (a) 
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to  the  gas  heating  term  Ah  are  ion  heating,  elastic  collisions 
and  rotational  excitation,  and  30%  of  the  power  deposited 
in  electronic  excitation,  i.e.  Ah  =  (Ai-r+A)  +  Aons 
(vibrational  to  translational  energy  transfer  is  not  included). 

Figure  7  compares  the  pressure  wave  calculated  at  time 
t  =  8 /as,  for  four  different  assumptions  on  the  gas  heating 
mechanisms:  no  gas  heating  i.e.  the  only  effect  of  the  plasma 
on  the  flow  is  the  EHD  force  (Ah  =  0),  ion  contribution 
only  (Pth  =  Aons).  contribution  of  ions,  elastic  collisions 
and  rotational  excitation,  and  quenching  of  the  excited  states 
(Ah  =  (Pei-R  +  Pe)  +  Aons),  and  same  plus  contribution 
of  vibrational  to  translational  (V-T)  energy  transfer  (Ah  = 
(Pei— r  +  Pe  +  At)  +  Aons),  with  tvt  =  1/zs.  We  see  in 
this  figure  that  the  amplitude  and  the  shape  of  the  micro 
shock  wave  are  significantly  different  when  one  or  the  other 
of  these  four  assumptions  is  made.  The  micro  shock  wave 
is  practically  absent  in  the  first  case  (figure  7(a))  where  gas 
heating  is  not  included.  The  small  effect  that  can  be  observed 


Figure  8.  M  aximum  gas  temperature  as  a  function  of  time  for 
different  assumptions  on  gas  heating. 

in  figure  7(a)  is  due  to  the  EHD  force.  The  amplitude  of 
the  pressure  wave  increases  by  a  factor  of  about  4  when  the 
contribution  of  quenching  of  the  excited  states  is  included,  and 
by  another  factor  of  2  when  V-T  energy  transfer  is  taken  into 
account. 

The  cylindrical  part  of  the  pressure  perturbation  that  can 
be  observed  in  figures  6  and  7  is  due  to  the  fact  that  a  large  part 
of  the  energy  i  s  rel  eased  at  the  ti  p  of  the  top  el  ectrode  duri  ng  the 
second  pulse,  i.e.  in  the  cathode  region  of  the  discharge  where 
large  ion  heating  and  electron  excitation  take  place.  The  more 
complex  shape  around  the  overall  cylindrical  shape  is  due  to 
the  fact  that  a  large  part  of  the  energy  is  also  released  during  the 
filament  propagation  along  the  surface  (in  the  streamer  head 
during  the  first  pulse,  and  in  the  plasma  channel).  This  is  the 
reason  why  the  pressure  perturbation  of  figures  6  and  7  has  a 
large  component  parallel  to  the  surface. 

The  propagation  velocity  of  the  shock  wave  can  be 
deduced  from  the  pressure  distribution  at  different  times.  We 
find  that,  in  the  conditions  of  figure  6,  the  propagation  velocity 
of  the  micro  shock  wave  is  450  m  s-1  during  the  first  100  ns, 
and  quickly  decreases  to  about  350 ms-1,  close  to  the  sound 
velocity.  This  is  in  good  agreement  with  the  experimental 
results  of  Stari kovski i  et  al  [18]  except  that  the  decay  of  the 
velocity  occurs  at  later  times  (on  the  order  of  1  /xs)  in  the 
experiments. 

The  micro  shock  wave  generation  is  due  to  the  very 
fast  energy  release  into  gas  heating.  We  performed  similar 
calculations  for  si  nusoidal  voltages  i  n  the  1- 10  kH  z  range.  T  he 
streamers  that  form  under  these  conditions  are  less  energetic 
(because  of  the  slower  voltage  rise,  the  overvoltage  is  not  as 
important  and  the  streamer  forms  at  a  lower  voltage,  closer 
to  the  breakdown  value)  and  the  current  pulse  duration  is 
longer.  A  smaller  energy  is  deposited  on  a  longer  timescale 
and  therefore  no  significant  shock  wave  is  observed  in  the 
simulation  under  these  conditions. 

Figure  8  shows  the  time  evolution  of  the  calculated 
maximum  gas  temperature  for  different  assumptions  on  the 
heating  term.  We  see  that  the  maximum  gas  temperature 
can  be  as  large  as  1000  K  at  the  end  of  the  discharge  pulses, 
under  these  conditions.  The  gas  temperature  comes  back 
to  ambient  temperature  on  a  microsecond  timescale.  As 
expected,  becauseof  the  time  del  ay  in  theV-T  energy  transfer, 
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Figure  9.  D  istribution  of  the  gas  temperature  after  the  end  of  the 
discharge  (timer  =  40  ns),  under  the  conditions  of  figures  3-5. 

the  maximum  temperature  is  not  (but  the  temperature  decay  is) 
affected  by  theV-T  transfer. 

The  space  distribution  of  the  gas  temperature  at  time 
t  =  40  ns,  just  after  the  discharge  pulse  is  shown  in  figure  9. 
We  clearly  see  a  hot  spot  in  the  region  close  to  the  tip  of  the 
top  electrode,  where  a  cathode  sheath  forms  during  the  second 
current  pulse.  The  power  density  deposited  in  gas  heating 
in  this  cathode  sheath  region  is  quite  large  and  leads  to  fast 
heating.  The  energy  deposited  during  the  first  current  pulse  is 
also  large  (even  larger  than  during  the  second  current  pulse, 
see  below),  but  is  distributed  over  a  larger  volume,  during  the 
streamer  propagation  along  the  surface.  This  is  the  reason  why 
the  gas  temperature  above  the  electrode  above  the  dielectric 
layer  (x  >  3.5  cm)  is  smaller,  on  the  order  of  400-500  K,  but 
the  temperature  increase  takes  place  in  a  larger  volume. 

Calculations  of  the  maximum  gas  temperature  have  also 
been  performed  for  different  values  of  the  minimum  cell  size. 
The  results  show  that  the  calculated  maximum  temperature 
can  be  significantly  larger  than  1000  K  when  the  minimum 
cell  size  is  5  urn  instead  of  19 ^m.  This  difference  is  due  to 
the  poor  resolution  of  the  cathode  sheath  when  the  mesh  is  too 
large.  The  local  power  density  is  significantly  larger  for  the 
smallest  mesh,  but  the  total  energy  deposition  is  only  10-20% 
larger  when  the  minimum  cell  size  is  changed  from  19  to  5  fim. 
The  calculated  pressure  wave  amplitude  at  t  =  8 /xs  is  about 
twice  as  large  for  the  smallest  mesh.  Note  that  it  is  difficult 
to  get  a  very  accurate  estimation  of  the  maximum  local  power 
deposition  not  only  because  this  needs  a  very  fine  resolution 
of  the  cathode  or  streamer  sheath  but  also  because  the  physical 
model  is  limited  by  the  local  field  approximation  (i.e.  the 
electron  energy  deposition  by  electrons  in  the  sheath  actually 
takes  place  over  a  larger  region  because  of  non-local  effects; 
these  effects  tend  to  lower  the  maximum  power  deposition). 

The  energy  deposition  through  different  channels  is 
represented  in  figure  10  as  a  function  of  time.  M  ore  energy  is 
deposited  during  the  first  pulse,  since,  as  we  mentioned  earlier, 
the  second  pul se  transfers  a  smal I er  charge  than  the  first  one  in 
th i  s  si  m  u I  ati  o n .  W e  see  that  the  en ergy  rel  ease  f  ro m  el  ectron i  c 
excitation  is  significantly  larger  than  the  energy  deposited  by 
ions.  The  largest  contribution  to  gas  heating  is  stored  in 
vibrational  excitation  but  is  released  on  a  longer  timescale. 

Figure  10  also  shows  that  the  total  energy  deposited  by 
the  discharge  is  on  the  order  of  0.5  mj  cm-1.  This  is  close  to 
the  experimental  value  given  by  Starikovskii  etal  [16,18]  for 
a  similar  voltage  pulse. 


Figure  10.  Calculated  energy  deposition  as  a  function  of  time: 
contribution  of  ions  to  gas  heating,  contribution  of  electron  elastic 
and  rotational  excitation  and  electronic  excitation  (30%  supposed  to 
be  converted  into  heating),  energy  stored  in  the  vibration  molecules 
and  total  energy  deposited  by  the  discharge  in  the  gas. 

5.  Conclusion 

We  have  developed  a  model  to  describe  plasma-flow 
interaction  under  conditions  of  fast  energy  deposition  by 
a  surface  discharge  plasma  (nanosecond  voltage  pulse  in  a 
SDBD),  leading  to  shock  wave  formation  and  propagation. 
The  model  solves  fluid  discharge  equations  coupled  with  the 
compressible  form  of  Navier- Stokes  equations. 

Thenumerical  model  isbased  on  an  original  asynchronous 
time  integration  method  combined  with  an  adaptative  mesh 
refinementtechnique.  This  asynchronous  time  integration  uses 
a  local  CFL  constraint,  so  less  computation  time  is  needed  in 
less  active  regions  where  the  particle  velocity  is  smaller  or 
the  grid  spacing  is  larger.  Because  of  the  adaptative  mesh 
refinement,  theAAM  R  method  also  allows  to  cope  with  large 
simulation  domains.  The  method  has  proven  to  be  very 
efficient  (the  speed  of  the  simulation  can  be  increased  by  a 
factor  of  100  with  respect  to  standard  methods),  provided 
that  care  is  taken  in  the  coding  and  in  the  algorithms  for 
asynchronous  timeinteg  ration  and  adaptative  mesh  refinement 
(the  details  are  not  described  here). 

The  simulations  have  been  performed  for  a  14  kV  voltage 
pulse  of  rise  time  7  ns,  duration  (plateau)  13  ns  and  decay  15  ns, 
and  for  a  discharge  over  a  300  fim  thick  dielectric  layer  and 
with  5  mm  long  electrodes  in  a  standard  SDBD  geometry.  The 
gas  dynamics  induced  by  gas  heating  is  followed  for  several 
microseconds  after  the  end  of  the  voltage  pulse. 

The  following  conclusions  can  be  drawn: 

•  A  first  current  pulse  takes  place  during  the  voltage  rise, 
when  the  top  electrode  plays  the  role  of  an  anode.  During 
this  pulsea  streamer  dischargeforms  and  spreads  along  the 
dielectric  surface.  Thestreamer  head  propagates  along  the 
dielectric  surface  at  a  velocity  of  around  5  x  107  cm  s-1. 
The  density  in  the  streamer  channel  is  on  the  order  of 
1020  m3.  During  the  fast  decay  of  the  applied  voltage,  the 
potential  difference  between  the  charged  dielectric  surface 
and  the  top  electrode  quickly  increases  to  values  over  the 
breakdown  voltage  and  a  second  current  pulse  is  initiated. 
During  this  pulse  the  exposed  electrode  plays  the  role  of 
a  cathode  and  a  cathode  sheath  forms  over  it. 
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•  The  fast  energy  deposition  under  the  considered 
conditions  of  a  nanosecond  voltage  pulse  leads  to  a  fast 
increase  in  the  gas  temperature  in  the  discharge  volume, 
with  a  very  localized  maximum  of  about  1000  K  next  to 
the  ti  p  of  the  exposed  el  ectrode,  at  the  end  of  the  di  scharge 
pulse.  The  location  of  this  maximum  coincides  with 
the  cathode  sheath  region  of  the  second  discharge  pulse. 
This  high  gas  temperature  is  associated  with  a  very  large 
pressure  increase  in  a  small  volume. 

•  The  local  pressure  increase  due  to  fast  gas  heating  leads  to 
the  generation  and  propagation  of  a  micro  shock  wave  with 
a  pressure  increase  of  a  few  1000  Pa  when  the  size  of  the 
perturbation  is  in  the  millimetre  range.  The  magnitude  of 
the  pressure  increase  depends  on  the  detailed  mechanisms 
of  energy  transfer  from  electron  to  the  gas  molecules.  A 
small  effect  is  seen  when  only  instantaneous  ion  heating  is 
considered.  This  effect  is  multiplied  by  4  when  30%  of  the 
energy  deposited  in  electronic  excitation  by  electrons  is 
supposed  to  be  i  nstantaneously  converted  i  nto  gas  heati  ng, 
and  by  a  factor  of  8  when  the  contribution  of  V-T  energy 
transfer  is  included,  with  a  V-T  relaxation  time  of  1  /xs. 

Although  the  gas  heating  model  used  in  this  paper  is 
approximate,  the  presented  results  are  consistent  with  the 
experiments  of  Starikovskii  et  al  if  a  large  part  of  the  energy 
absorbed  by  the  electrons  is  supposed  to  be  quickly  released 
into  gas  heating.  The  numerical  results  confirm  that  in  contrast 
with  the  sinusoidal  regime  of  DBDs  where  momentum  gained 
by  ions  from  the  field  is  transferred  to  the  gas  (ion  wind 
effect),  the  aerodynamics  effects  in  a  nanosecond  discharge 
are  associated  with  the  fast  gas  heating  in  the  boundary  layer. 
Whether  the  observed  effects  on  the  flow  are  due  to  the  gas 
dynamics  induced  by  the  nanosecond  discharge  or  to  changes 
in  the  viscosity  due  to  gas  heating  will  be  the  subject  of  further 
investigations. 
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